JETSCAPE/external_packages/fjcore.hh

0001 // fjcore -- extracted from FastJet v3.2.1 (http://fastjet.fr)
0002 //
0003 // fjcore constitutes a digest of the main FastJet functionality.
0004 // The files fjcore.hh and fjcore.cc are meant to provide easy access to these
0005 // core functions, in the form of single files and without the need of a full
0006 // FastJet installation:
0007 //
0008 //     g++ main.cc fjcore.cc
0009 //
0010 // with main.cc including fjcore.hh.
0011 //
0012 // A fortran interface, fjcorefortran.cc, is also provided. See the example
0013 // and the Makefile for instructions.
0014 //
0015 // The results are expected to be identical to those obtained by linking to
0016 // the full FastJet distribution.
0017 //
0018 // NOTE THAT, IN ORDER TO MAKE IT POSSIBLE FOR FJCORE AND THE FULL FASTJET
0019 // TO COEXIST, THE FORMER USES THE "fjcore" NAMESPACE INSTEAD OF "fastjet".
0020 //
0021 // In particular, fjcore provides:
0022 //
0023 //   - access to all native pp and ee algorithms, kt, anti-kt, C/A.
0024 //     For C/A, the NlnN method is available, while anti-kt and kt
0025 //     are limited to the N^2 one (still the fastest for N < 100k particles)
0026 //   - access to selectors, for implementing cuts and selections
0027 //   - access to all functionalities related to pseudojets (e.g. a jet's
0028 //     structure or user-defined information)
0029 //
0030 // Instead, it does NOT provide:
0031 //
0032 //   - jet areas functionality
0033 //   - background estimation
0034 //   - access to other algorithms via plugins
0035 //   - interface to CGAL
0036 //   - fastjet tools, e.g. filters, taggers
0037 //
0038 // If these functionalities are needed, the full FastJet installation must be
0039 // used. The code will be fully compatible, with the sole replacement of the
0040 // header files and of the fjcore namespace with the fastjet one.
0041 //
0042 // fjcore.hh and fjcore.cc are not meant to be human-readable.
0043 // For documentation, see the full FastJet manual and doxygen at http://fastjet.fr
0044 //
0045 // Like FastJet, fjcore is released under the terms of the GNU General Public
0046 // License version 2 (GPLv2). If you use this code as part of work towards a
0047 // scientific publication, whether directly or contained within another program
0048 // (e.g. Delphes, MadGraph, SpartyJet, Rivet, LHC collaboration software frameworks,
0049 // etc.), you should include a citation to
0050 //
0051 //   EPJC72(2012)1896 [arXiv:1111.6097] (FastJet User Manual)
0052 //   and, optionally, Phys.Lett.B641 (2006) 57 [arXiv:hep-ph/0512210]
0053 //
0054 // Copyright (c) 2005-2016, Matteo Cacciari, Gavin P. Salam and Gregory Soyez
0055 //
0056 //----------------------------------------------------------------------
0057 // This file is part of FastJet.
0058 //
0059 //  FastJet is free software; you can redistribute it and/or modify
0060 //  it under the terms of the GNU General Public License as published by
0061 //  the Free Software Foundation; either version 2 of the License, or
0062 //  (at your option) any later version.
0063 //
0064 //  The algorithms that underlie FastJet have required considerable
0065 //  development and are described in hep-ph/0512210. If you use
0066 //  FastJet as part of work towards a scientific publication, please
0067 //  include a citation to the FastJet paper.
0068 //
0069 //  FastJet is distributed in the hope that it will be useful,
0070 //  but WITHOUT ANY WARRANTY; without even the implied warranty of
0071 //  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
0072 //  GNU General Public License for more details.
0073 //
0074 //  You should have received a copy of the GNU General Public License
0075 //  along with FastJet. If not, see <http://www.gnu.org/licenses/>.
0076 //----------------------------------------------------------------------
0077 //
0078 #ifndef __FJCORE_HH__
0079 #define __FJCORE_HH__
0080 #define __FJCORE__   // remove all the non-core code (a safekeeper)
0081 #define __FJCORE_DROP_CGAL    // disable CGAL support
0082 #ifndef _INCLUDE_FJCORE_CONFIG_AUTO_H
0083 #define _INCLUDE_FJCORE_CONFIG_AUTO_H 1
0084 #ifndef FJCORE_HAVE_DEMANGLING_SUPPORT
0085 #endif
0086 #ifndef FJCORE_HAVE_DLFCN_H
0087 # define FJCORE_HAVE_DLFCN_H  1
0088 #endif
0089 #ifndef FJCORE_HAVE_EXECINFO_H
0090 #endif
0091 #ifndef FJCORE_HAVE_GNUCXX_DEPRECATED
0092 #endif
0093 #ifndef FJCORE_HAVE_INTTYPES_H
0094 # define FJCORE_HAVE_INTTYPES_H  1
0095 #endif
0096 #ifndef FJCORE_HAVE_LIBM
0097 # define FJCORE_HAVE_LIBM  1
0098 #endif
0099 #ifndef FJCORE_HAVE_MEMORY_H
0100 # define FJCORE_HAVE_MEMORY_H  1
0101 #endif
0102 #ifndef FJCORE_HAVE_STDINT_H
0103 # define FJCORE_HAVE_STDINT_H  1
0104 #endif
0105 #ifndef FJCORE_HAVE_STDLIB_H
0106 # define FJCORE_HAVE_STDLIB_H  1
0107 #endif
0108 #ifndef FJCORE_HAVE_STRINGS_H
0109 # define FJCORE_HAVE_STRINGS_H  1
0110 #endif
0111 #ifndef FJCORE_HAVE_STRING_H
0112 # define FJCORE_HAVE_STRING_H  1
0113 #endif
0114 #ifndef FJCORE_HAVE_SYS_STAT_H
0115 # define FJCORE_HAVE_SYS_STAT_H  1
0116 #endif
0117 #ifndef FJCORE_HAVE_SYS_TYPES_H
0118 # define FJCORE_HAVE_SYS_TYPES_H  1
0119 #endif
0120 #ifndef FJCORE_HAVE_UNISTD_H
0121 # define FJCORE_HAVE_UNISTD_H  1
0122 #endif
0123 #ifndef FJCORE_LT_OBJDIR
0124 # define FJCORE_LT_OBJDIR  ".libs/"
0125 #endif
0126 #ifndef FJCORE_PACKAGE
0127 # define FJCORE_PACKAGE  "fastjet"
0128 #endif
0129 #ifndef FJCORE_PACKAGE_BUGREPORT
0130 # define FJCORE_PACKAGE_BUGREPORT  ""
0131 #endif
0132 #ifndef FJCORE_PACKAGE_NAME
0133 # define FJCORE_PACKAGE_NAME  "FastJet"
0134 #endif
0135 #ifndef FJCORE_PACKAGE_STRING
0136 # define FJCORE_PACKAGE_STRING  "FastJet 3.2.1"
0137 #endif
0138 #ifndef FJCORE_PACKAGE_TARNAME
0139 # define FJCORE_PACKAGE_TARNAME  "fastjet"
0140 #endif
0141 #ifndef FJCORE_PACKAGE_URL
0142 # define FJCORE_PACKAGE_URL  ""
0143 #endif
0144 #ifndef FJCORE_PACKAGE_VERSION
0145 # define FJCORE_PACKAGE_VERSION  "3.2.1"
0146 #endif
0147 #ifndef FJCORE_STDC_HEADERS
0148 # define FJCORE_STDC_HEADERS  1
0149 #endif
0150 #ifndef FJCORE_VERSION
0151 # define FJCORE_VERSION  "3.2.1"
0152 #endif
0153 #ifndef FJCORE_VERSION_MAJOR
0154 # define FJCORE_VERSION_MAJOR  3
0155 #endif
0156 #ifndef FJCORE_VERSION_MINOR
0157 # define FJCORE_VERSION_MINOR  2
0158 #endif
0159 #ifndef FJCORE_VERSION_NUMBER
0160 # define FJCORE_VERSION_NUMBER  30201
0161 #endif
0162 #ifndef FJCORE_VERSION_PATCHLEVEL
0163 # define FJCORE_VERSION_PATCHLEVEL  1
0164 #endif
0165 #endif
0166 #ifndef __FJCORE_CONFIG_H__
0167 #define __FJCORE_CONFIG_H__
0168 #endif // __FJCORE_CONFIG_H__
0169 #ifndef __FJCORE_FASTJET_BASE_HH__
0170 #define __FJCORE_FASTJET_BASE_HH__
0171 #define FJCORE_BEGIN_NAMESPACE namespace fjcore {
0172 #define FJCORE_END_NAMESPACE   }
0173 #ifdef FJCORE_HAVE_OVERRIDE
0174 # define FJCORE_OVERRIDE  override
0175 #else
0176 # define FJCORE_OVERRIDE
0177 #endif
0178 #endif // __FJCORE_FASTJET_BASE_HH__
0179 #ifndef __FJCORE_NUMCONSTS__
0180 #define __FJCORE_NUMCONSTS__
0181 FJCORE_BEGIN_NAMESPACE      // defined in fastjet/internal/base.hh
0182 const double pi = 3.141592653589793238462643383279502884197;
0183 const double twopi = 6.283185307179586476925286766559005768394;
0184 const double pisq  = 9.869604401089358618834490999876151135314;
0185 const double zeta2 = 1.644934066848226436472415166646025189219;
0186 const double zeta3 = 1.202056903159594285399738161511449990765;
0187 const double eulergamma = 0.577215664901532860606512090082402431042;
0188 const double ln2   = 0.693147180559945309417232121458176568076;
0189 FJCORE_END_NAMESPACE
0190 #endif // __FJCORE_NUMCONSTS__
0191 #ifndef __FJCORE_INTERNAL_IS_BASE_HH__
0192 #define __FJCORE_INTERNAL_IS_BASE_HH__
0193 FJCORE_BEGIN_NAMESPACE
0194 template<typename T, T _t>
0195 struct integral_type{
0196   static const T value = _t;         ///< the value (only member carrying info)
0197   typedef T value_type;          ///< a typedef for the type T
0198   typedef integral_type<T,_t> type;  ///< a typedef for the whole structure
0199 };
0200 template<typename T, T _t>
0201 const T integral_type<T, _t>::value;
0202 typedef integral_type<bool, true>  true_type;  ///< the bool 'true'  value promoted to a type
0203 typedef integral_type<bool, false> false_type; ///< the bool 'false' value promoted to a type
0204 typedef char (&__yes_type)[1]; //< the yes type
0205 typedef char (&__no_type) [2]; //< the no type
0206 template<typename B, typename D>
0207 struct __inheritance_helper{
0208 #if !((_MSC_VER !=0 ) && (_MSC_VER == 1310))   // MSVC 7.1
0209   template <typename T>
0210   static __yes_type check_sig(D const volatile *, T);
0211 #else
0212   static __yes_type check_sig(D const volatile *, long);
0213 #endif
0214   static __no_type  check_sig(B const volatile *, int);
0215 };
0216 template<typename B, typename D>
0217 struct IsBaseAndDerived{
0218 #if ((_MSC_FULL_VER != 0) && (_MSC_FULL_VER >= 140050000))
0219 #pragma warning(push)
0220 #pragma warning(disable:6334)
0221 #endif
0222   struct Host{
0223 #if !((_MSC_VER !=0 ) && (_MSC_VER == 1310))
0224     operator B const volatile *() const;
0225 #else
0226     operator B const volatile * const&() const;
0227 #endif
0228     operator D const volatile *();
0229   };
0230   static const bool value = ((sizeof(B)!=0) &&
0231                  (sizeof(D)!=0) &&
0232                  (sizeof(__inheritance_helper<B,D>::check_sig(Host(), 0)) == sizeof(__yes_type)));
0233 #if ((_MSC_FULL_VER != 0) && (_MSC_FULL_VER >= 140050000))
0234 #pragma warning(pop)
0235 #endif
0236 };
0237 template<class B, class D>
0238 B* cast_if_derived(D* d){
0239   return IsBaseAndDerived<B,D>::value ? (B*)(d) : 0;
0240 }
0241 FJCORE_END_NAMESPACE
0242 #endif  // __IS_BASE_OF_HH__
0243 #ifndef __FJCORE_FJCORE_DEPRECATED_HH__
0244 #define __FJCORE_FJCORE_DEPRECATED_HH__
0245 #if defined(FJCORE_HAVE_CXX14_DEPRECATED)
0246 # define FJCORE_DEPRECATED               [[deprecated]]
0247 # define FJCORE_DEPRECATED_MSG(message)  [[deprecated(message)]]
0248 #elif defined(FJCORE_HAVE_GNUCXX_DEPRECATED)
0249 # define FJCORE_DEPRECATED               __attribute__((__deprecated__))
0250 # define FJCORE_DEPRECATED_MSG(message)  __attribute__((__deprecated__))
0251 #else
0252 # define FJCORE_DEPRECATED
0253 # define FJCORE_DEPRECATED_MSG(message)
0254 #endif
0255 #endif // __FJCORE_FJCORE_DEPRECATED_HH__
0256 #ifndef __FJCORE_SHARED_PTR_HH__
0257 #define __FJCORE_SHARED_PTR_HH__
0258 #include <cstdlib>  // for NULL!!!
0259 #ifdef __FJCORE_USETR1SHAREDPTR
0260 #include <tr1/memory>
0261 #endif // __FJCORE_USETR1SHAREDPTR
0262 FJCORE_BEGIN_NAMESPACE      // defined in fastjet/internal/base.hh
0263 #ifdef __FJCORE_USETR1SHAREDPTR
0264 template<class T>
0265 class SharedPtr : public std::tr1::shared_ptr<T> {
0266 public:
0267   SharedPtr() : std::tr1::shared_ptr<T>() {}
0268   SharedPtr(T * t) : std::tr1::shared_ptr<T>(t) {}
0269   SharedPtr(const SharedPtr<T> & t) : std::tr1::shared_ptr<T>(t) {}
0270   #ifdef FJCORE_HAVE_EXPLICIT_FOR_OPERATORS
0271   explicit
0272   #endif
0273   inline operator bool() const {return (this->get()!=NULL);}
0274   T* operator ()() const{
0275     return this->get(); // automatically returns NULL when out-of-scope
0276   }
0277 };
0278 #else // __FJCORE_USETR1SHAREDPTR
0279 template<class T>
0280 class SharedPtr{
0281 public:
0282   class __SharedCountingPtr;
0283   SharedPtr() : _ptr(NULL){}
0284   template<class Y> explicit SharedPtr(Y* ptr){
0285     _ptr = new __SharedCountingPtr(ptr);
0286   }
0287   SharedPtr(SharedPtr const & share) : _ptr(share._get_container()){
0288     if (_ptr!=NULL) ++(*_ptr);
0289   }
0290   ~SharedPtr(){
0291     if (_ptr==NULL) return;
0292     _decrease_count();
0293   }
0294   void reset(){
0295     SharedPtr().swap(*this);
0296   }
0297   template<class Y> void reset(Y * ptr){
0298     SharedPtr(ptr).swap(*this);
0299   }
0300   template<class Y> void reset(SharedPtr<Y> const & share){
0301     if (_ptr!=NULL){
0302       if (_ptr == share._get_container()) return;
0303       _decrease_count();
0304     }
0305     _ptr = share._get_container();  // Note: automatically set it to NULL if share is empty
0306     if (_ptr!=NULL) ++(*_ptr);
0307   }
0308   SharedPtr& operator=(SharedPtr const & share){
0309     reset(share);
0310     return *this;
0311   }
0312   template<class Y> SharedPtr& operator=(SharedPtr<Y> const & share){
0313     reset(share);
0314     return *this;
0315   }
0316   FJCORE_DEPRECATED_MSG("Use SharedPtr<T>::get() instead")
0317   T* operator ()() const{
0318     if (_ptr==NULL) return NULL;
0319     return _ptr->get(); // automatically returns NULL when out-of-scope
0320   }
0321   inline T& operator*() const{
0322     return *(_ptr->get());
0323   }
0324   inline T* operator->() const{
0325     if (_ptr==NULL) return NULL;
0326     return _ptr->get();
0327   }
0328   inline T* get() const{
0329     if (_ptr==NULL) return NULL;
0330     return _ptr->get();
0331   }
0332   inline bool unique() const{
0333     return (use_count()==1);
0334   }
0335   inline long use_count() const{
0336     if (_ptr==NULL) return 0;
0337     return _ptr->use_count(); // automatically returns NULL when out-of-scope
0338   }
0339   #ifdef FJCORE_HAVE_EXPLICIT_FOR_OPERATORS
0340   explicit
0341   #endif
0342   inline operator bool() const{
0343     return (get()!=NULL);
0344   }
0345   inline void swap(SharedPtr & share){
0346     __SharedCountingPtr* share_container = share._ptr;
0347     share._ptr = _ptr;
0348     _ptr = share_container;
0349   }
0350   void set_count(const long & count){
0351     if (_ptr==NULL) return;
0352     _ptr->set_count(count);
0353   }
0354   class __SharedCountingPtr{
0355   public:
0356     __SharedCountingPtr() : _ptr(NULL), _count(0){}
0357     template<class Y> explicit __SharedCountingPtr(Y* ptr) : _ptr(ptr), _count(1){}
0358     ~__SharedCountingPtr(){
0359       if (_ptr!=NULL){ delete _ptr;}
0360     }
0361     inline T* get() const {return _ptr;}
0362     inline long use_count() const {return _count;}
0363     inline long operator++(){return ++_count;}
0364     inline long operator--(){return --_count;}
0365     inline long operator++(int){return _count++;}
0366     inline long operator--(int){return _count--;}
0367     void set_count(const long & count){
0368       _count = count;
0369     }
0370   private:
0371     T *_ptr;              ///< the pointer we're counting the references to
0372     long _count;  ///< the number of references
0373   };
0374 private:
0375   inline __SharedCountingPtr* _get_container() const{
0376     return _ptr;
0377   }
0378   void _decrease_count(){
0379     (*_ptr)--;
0380     if (_ptr->use_count()==0)
0381       delete _ptr; // that automatically deletes the object itself
0382   }
0383   __SharedCountingPtr *_ptr;
0384 };
0385 template<class T,class U>
0386 inline bool operator==(SharedPtr<T> const & t, SharedPtr<U> const & u){
0387   return t.get() == u.get();
0388 }
0389 template<class T,class U>
0390 inline bool operator!=(SharedPtr<T> const & t, SharedPtr<U> const & u){
0391   return t.get() != u.get();
0392 }
0393 template<class T,class U>
0394 inline bool operator<(SharedPtr<T> const & t, SharedPtr<U> const & u){
0395   return t.get() < u.get();
0396 }
0397 template<class T>
0398 inline void swap(SharedPtr<T> & a, SharedPtr<T> & b){
0399   return a.swap(b);
0400 }
0401 template<class T>
0402 inline T* get_pointer(SharedPtr<T> const & t){
0403   return t.get();
0404 }
0405 #endif // __FJCORE_USETR1SHAREDPTR
0406 FJCORE_END_NAMESPACE      // defined in fastjet/internal/base.hh
0407 #endif   // __FJCORE_SHARED_PTR_HH__
0408 #ifndef __FJCORE_LIMITEDWARNING_HH__
0409 #define __FJCORE_LIMITEDWARNING_HH__
0410 #include <iostream>
0411 #include <string>
0412 #include <list>
0413 FJCORE_BEGIN_NAMESPACE
0414 class LimitedWarning {
0415 public:
0416   LimitedWarning() : _max_warn(_max_warn_default), _n_warn_so_far(0), _this_warning_summary(0) {}
0417   LimitedWarning(int max_warn_in) : _max_warn(max_warn_in), _n_warn_so_far(0), _this_warning_summary(0) {}
0418   void warn(const char * warning) {warn(warning, _default_ostr);}
0419   void warn(const std::string & warning) {warn(warning.c_str(), _default_ostr);}
0420   void warn(const char * warning, std::ostream * ostr);
0421   void warn(const std::string & warning, std::ostream * ostr) {warn(warning.c_str(), ostr);}
0422   static void set_default_stream(std::ostream * ostr) {
0423     _default_ostr = ostr;
0424   }
0425   static void set_default_max_warn(int max_warn) {
0426     _max_warn_default = max_warn;
0427   }
0428   int max_warn() const {return _max_warn;}
0429   int n_warn_so_far() const {return _n_warn_so_far;}
0430   static std::string summary();
0431 private:
0432   int _max_warn, _n_warn_so_far;
0433   static int _max_warn_default;
0434   static std::ostream * _default_ostr;
0435   typedef std::pair<std::string, unsigned int> Summary;
0436   static std::list< Summary > _global_warnings_summary;
0437   Summary * _this_warning_summary;
0438 };
0439 FJCORE_END_NAMESPACE
0440 #endif // __FJCORE_LIMITEDWARNING_HH__
0441 #ifndef __FJCORE_ERROR_HH__
0442 #define __FJCORE_ERROR_HH__
0443 #include<iostream>
0444 #include<string>
0445 #if (!defined(FJCORE_HAVE_EXECINFO_H)) || defined(__FJCORE__)
0446 #endif
0447 FJCORE_BEGIN_NAMESPACE      // defined in fastjet/internal/base.hh
0448 class Error {
0449 public:
0450   Error() {}
0451   Error(const std::string & message);
0452   virtual ~Error() {}
0453   std::string message() const {return _message;}
0454   static void set_print_errors(bool print_errors) {_print_errors = print_errors;}
0455   static void set_print_backtrace(bool enabled);
0456   static void set_default_stream(std::ostream * ostr) {
0457     _default_ostr = ostr;
0458   }
0459 private:
0460   std::string _message;                ///< error message
0461   static bool _print_errors;           ///< do we print anything?
0462   static bool _print_backtrace;        ///< do we print the backtrace?
0463   static std::ostream * _default_ostr; ///< the output stream (cerr if not set)
0464 #if (!defined(FJCORE_HAVE_EXECINFO_H)) || defined(__FJCORE__)
0465   static LimitedWarning _execinfo_undefined;
0466 #endif
0467 };
0468 class InternalError : public Error{
0469 public:
0470   InternalError(const std::string & message_in) : Error(std::string("*** CRITICAL INTERNAL FASTJET ERROR *** CONTACT THE AUTHORS *** ") + message_in){ }
0471 };
0472 FJCORE_END_NAMESPACE
0473 #endif // __FJCORE_ERROR_HH__
0474 #ifndef __FJCORE_PSEUDOJET_STRUCTURE_BASE_HH__
0475 #define __FJCORE_PSEUDOJET_STRUCTURE_BASE_HH__
0476 #include <vector>
0477 #include <string>
0478 FJCORE_BEGIN_NAMESPACE      // defined in fastjet/internal/base.hh
0479 class PseudoJet;
0480 class ClusterSequence;
0481 class PseudoJetStructureBase{
0482 public:
0483   PseudoJetStructureBase(){};
0484   virtual ~PseudoJetStructureBase(){};
0485   virtual std::string description() const{ return "PseudoJet with an unknown structure"; }
0486   virtual bool has_associated_cluster_sequence() const { return false;}
0487   virtual const ClusterSequence* associated_cluster_sequence() const;
0488   virtual bool has_valid_cluster_sequence() const {return false;}
0489   virtual const ClusterSequence * validated_cs() const;
0490   virtual bool has_partner(const PseudoJet &reference, PseudoJet &partner) const;
0491   virtual bool has_child(const PseudoJet &reference, PseudoJet &child) const;
0492   virtual bool has_parents(const PseudoJet &reference, PseudoJet &parent1, PseudoJet &parent2) const;
0493   virtual bool object_in_jet(const PseudoJet &reference, const PseudoJet &jet) const;
0494   virtual bool has_constituents() const {return false;}
0495   virtual std::vector<PseudoJet> constituents(const PseudoJet &reference) const;
0496   virtual bool has_exclusive_subjets() const {return false;}
0497   virtual std::vector<PseudoJet> exclusive_subjets(const PseudoJet &reference, const double & dcut) const;
0498   virtual int n_exclusive_subjets(const PseudoJet &reference, const double & dcut) const;
0499   virtual std::vector<PseudoJet> exclusive_subjets_up_to (const PseudoJet &reference, int nsub) const;
0500   virtual double exclusive_subdmerge(const PseudoJet &reference, int nsub) const;
0501   virtual double exclusive_subdmerge_max(const PseudoJet &reference, int nsub) const;
0502   virtual bool has_pieces(const PseudoJet & /* reference */) const {
0503     return false;}
0504   virtual std::vector<PseudoJet> pieces(const PseudoJet & /* reference */
0505                                         ) const;
0506 };
0507 FJCORE_END_NAMESPACE
0508 #endif  //  __FJCORE_PSEUDOJET_STRUCTURE_BASE_HH__
0509 #ifndef __FJCORE_PSEUDOJET_HH__
0510 #define __FJCORE_PSEUDOJET_HH__
0511 #include<valarray>
0512 #include<vector>
0513 #include<cassert>
0514 #include<cmath>
0515 #include<iostream>
0516 FJCORE_BEGIN_NAMESPACE      // defined in fastjet/internal/base.hh
0517 const double MaxRap = 1e5;
0518 const double pseudojet_invalid_phi = -100.0;
0519 const double pseudojet_invalid_rap = -1e200;
0520 class PseudoJet {
0521  public:
0522   PseudoJet() : _px(0), _py(0), _pz(0), _E(0) {_finish_init(); _reset_indices();}
0523   PseudoJet(const double px, const double py, const double pz, const double E);
0524   template <class L> PseudoJet(const L & some_four_vector);
0525   PseudoJet(bool /* dummy */) {}
0526   virtual ~PseudoJet(){};
0527   inline double E()   const {return _E;}
0528   inline double e()   const {return _E;} // like CLHEP
0529   inline double px()  const {return _px;}
0530   inline double py()  const {return _py;}
0531   inline double pz()  const {return _pz;}
0532   inline double phi() const {return phi_02pi();}
0533   inline double phi_std()  const {
0534     _ensure_valid_rap_phi();
0535     return _phi > pi ? _phi-twopi : _phi;}
0536   inline double phi_02pi() const {
0537     _ensure_valid_rap_phi();
0538     return _phi;
0539   }
0540   inline double rap() const {
0541     _ensure_valid_rap_phi();
0542     return _rap;
0543   }
0544   inline double rapidity() const {return rap();} // like CLHEP
0545   double pseudorapidity() const;
0546   double eta() const {return pseudorapidity();}
0547   inline double pt2() const {return _kt2;}
0548   inline double  pt() const {return sqrt(_kt2);}
0549   inline double perp2() const {return _kt2;}  // like CLHEP
0550   inline double  perp() const {return sqrt(_kt2);}    // like CLHEP
0551   inline double kt2() const {return _kt2;} // for bkwds compatibility
0552   inline double  m2() const {return (_E+_pz)*(_E-_pz)-_kt2;}
0553   inline double  m() const;
0554   inline double mperp2() const {return (_E+_pz)*(_E-_pz);}
0555   inline double mperp() const {return sqrt(std::abs(mperp2()));}
0556   inline double mt2() const {return (_E+_pz)*(_E-_pz);}
0557   inline double mt() const {return sqrt(std::abs(mperp2()));}
0558   inline double modp2() const {return _kt2+_pz*_pz;}
0559   inline double modp() const {return sqrt(_kt2+_pz*_pz);}
0560   inline double Et() const {return (_kt2==0) ? 0.0 : _E/sqrt(1.0+_pz*_pz/_kt2);}
0561   inline double Et2() const {return (_kt2==0) ? 0.0 : _E*_E/(1.0+_pz*_pz/_kt2);}
0562   double operator () (int i) const ;
0563   inline double operator [] (int i) const { return (*this)(i); }; // this too
0564   double kt_distance(const PseudoJet & other) const;
0565   double plain_distance(const PseudoJet & other) const;
0566   inline double squared_distance(const PseudoJet & other) const {
0567     return plain_distance(other);}
0568   inline double delta_R(const PseudoJet & other) const {
0569     return sqrt(squared_distance(other));
0570   }
0571   double delta_phi_to(const PseudoJet & other) const;
0572   inline double beam_distance() const {return _kt2;}
0573   std::valarray<double> four_mom() const;
0574   enum { X=0, Y=1, Z=2, T=3, NUM_COORDINATES=4, SIZE=NUM_COORDINATES };
0575   PseudoJet & boost(const PseudoJet & prest);
0576   PseudoJet & unboost(const PseudoJet & prest);
0577   void operator*=(double);
0578   void operator/=(double);
0579   void operator+=(const PseudoJet &);
0580   void operator-=(const PseudoJet &);
0581   inline void reset(double px, double py, double pz, double E);
0582   inline void reset(const PseudoJet & psjet) {
0583     (*this) = psjet;
0584   }
0585   template <class L> inline void reset(const L & some_four_vector) {
0586     const PseudoJet * pj = fjcore::cast_if_derived<const PseudoJet>(&some_four_vector);
0587     if (pj){
0588       (*this) = *pj;
0589     } else {
0590       reset(some_four_vector[0], some_four_vector[1],
0591         some_four_vector[2], some_four_vector[3]);
0592     }
0593   }
0594   inline void reset_PtYPhiM(double pt_in, double y_in, double phi_in, double m_in=0.0) {
0595     reset_momentum_PtYPhiM(pt_in, y_in, phi_in, m_in);
0596     _reset_indices();
0597   }
0598   inline void reset_momentum(double px, double py, double pz, double E);
0599   inline void reset_momentum(const PseudoJet & pj);
0600   void reset_momentum_PtYPhiM(double pt, double y, double phi, double m=0.0);
0601   template <class L> inline void reset_momentum(const L & some_four_vector) {
0602     reset_momentum(some_four_vector[0], some_four_vector[1],
0603            some_four_vector[2], some_four_vector[3]);
0604   }
0605   void set_cached_rap_phi(double rap, double phi);
0606   inline int user_index() const {return _user_index;}
0607   inline void set_user_index(const int index) {_user_index = index;}
0608   class UserInfoBase{
0609   public:
0610     UserInfoBase(){};
0611     virtual ~UserInfoBase(){};
0612   };
0613   class InexistentUserInfo : public Error {
0614   public:
0615     InexistentUserInfo();
0616   };
0617   void set_user_info(UserInfoBase * user_info_in) {
0618     _user_info.reset(user_info_in);
0619   }
0620   template<class L>
0621   const L & user_info() const{
0622     if (_user_info.get() == 0) throw InexistentUserInfo();
0623     return dynamic_cast<const L &>(* _user_info.get());
0624   }
0625   bool has_user_info() const{
0626     return _user_info.get();
0627   }
0628   template<class L>
0629   bool has_user_info() const{
0630     return _user_info.get() && dynamic_cast<const L *>(_user_info.get());
0631   }
0632   const UserInfoBase * user_info_ptr() const{
0633     return _user_info.get();
0634   }
0635   const SharedPtr<UserInfoBase> & user_info_shared_ptr() const{
0636     return _user_info;
0637   }
0638   SharedPtr<UserInfoBase> & user_info_shared_ptr(){
0639     return _user_info;
0640   }
0641   std::string description() const;
0642   bool has_associated_cluster_sequence() const;
0643   bool has_associated_cs() const {return has_associated_cluster_sequence();}
0644   bool has_valid_cluster_sequence() const;
0645   bool has_valid_cs() const {return has_valid_cluster_sequence();}
0646   const ClusterSequence* associated_cluster_sequence() const;
0647   const ClusterSequence* associated_cs() const {return associated_cluster_sequence();}
0648   inline const ClusterSequence * validated_cluster_sequence() const {
0649     return validated_cs();
0650   }
0651   const ClusterSequence * validated_cs() const;
0652   void set_structure_shared_ptr(const SharedPtr<PseudoJetStructureBase> &structure);
0653   bool has_structure() const;
0654   const PseudoJetStructureBase* structure_ptr() const;
0655   PseudoJetStructureBase* structure_non_const_ptr();
0656   const PseudoJetStructureBase* validated_structure_ptr() const;
0657   const SharedPtr<PseudoJetStructureBase> & structure_shared_ptr() const;
0658   template<typename StructureType>
0659   const StructureType & structure() const;
0660   template<typename TransformerType>
0661   bool has_structure_of() const;
0662   template<typename TransformerType>
0663   const typename TransformerType::StructureType & structure_of() const;
0664   virtual bool has_partner(PseudoJet &partner) const;
0665   virtual bool has_child(PseudoJet &child) const;
0666   virtual bool has_parents(PseudoJet &parent1, PseudoJet &parent2) const;
0667   virtual bool contains(const PseudoJet &constituent) const;
0668   virtual bool is_inside(const PseudoJet &jet) const;
0669   virtual bool has_constituents() const;
0670   virtual std::vector<PseudoJet> constituents() const;
0671   virtual bool has_exclusive_subjets() const;
0672   std::vector<PseudoJet> exclusive_subjets (const double dcut) const;
0673   int n_exclusive_subjets(const double dcut) const;
0674   std::vector<PseudoJet> exclusive_subjets (int nsub) const;
0675   std::vector<PseudoJet> exclusive_subjets_up_to (int nsub) const;
0676   double exclusive_subdmerge(int nsub) const;
0677   double exclusive_subdmerge_max(int nsub) const;
0678   virtual bool has_pieces() const;
0679   virtual std::vector<PseudoJet> pieces() const;
0680   inline int cluster_hist_index() const {return _cluster_hist_index;}
0681   inline void set_cluster_hist_index(const int index) {_cluster_hist_index = index;}
0682   inline int cluster_sequence_history_index() const {
0683     return cluster_hist_index();}
0684   inline void set_cluster_sequence_history_index(const int index) {
0685     set_cluster_hist_index(index);}
0686  protected:
0687   SharedPtr<PseudoJetStructureBase> _structure;
0688   SharedPtr<UserInfoBase> _user_info;
0689  private:
0690   double _px,_py,_pz,_E;
0691   mutable double _phi, _rap;
0692   double _kt2;
0693   int    _cluster_hist_index, _user_index;
0694   void _finish_init();
0695   void _reset_indices();
0696   inline void _ensure_valid_rap_phi() const {
0697     if (_phi == pseudojet_invalid_phi) _set_rap_phi();
0698   }
0699   void _set_rap_phi() const;
0700   friend PseudoJet operator*(double, const PseudoJet &);
0701 };
0702 PseudoJet operator+(const PseudoJet &, const PseudoJet &);
0703 PseudoJet operator-(const PseudoJet &, const PseudoJet &);
0704 PseudoJet operator*(double, const PseudoJet &);
0705 PseudoJet operator*(const PseudoJet &, double);
0706 PseudoJet operator/(const PseudoJet &, double);
0707 bool operator==(const PseudoJet &, const PseudoJet &);
0708 inline bool operator!=(const PseudoJet & a, const PseudoJet & b) {return !(a==b);}
0709 bool operator==(const PseudoJet & jet, const double val);
0710 inline bool operator==(const double val, const PseudoJet & jet) {return jet == val;}
0711 inline bool operator!=(const PseudoJet & a, const double val)  {return !(a==val);}
0712 inline bool operator!=( const double val, const PseudoJet & a) {return !(a==val);}
0713 inline double dot_product(const PseudoJet & a, const PseudoJet & b) {
0714   return a.E()*b.E() - a.px()*b.px() - a.py()*b.py() - a.pz()*b.pz();
0715 }
0716 bool have_same_momentum(const PseudoJet &, const PseudoJet &);
0717 PseudoJet PtYPhiM(double pt, double y, double phi, double m = 0.0);
0718 std::vector<PseudoJet> sorted_by_pt(const std::vector<PseudoJet> & jets);
0719 std::vector<PseudoJet> sorted_by_rapidity(const std::vector<PseudoJet> & jets);
0720 std::vector<PseudoJet> sorted_by_E(const std::vector<PseudoJet> & jets);
0721 std::vector<PseudoJet> sorted_by_pz(const std::vector<PseudoJet> & jets);
0722 void sort_indices(std::vector<int> & indices,
0723           const std::vector<double> & values);
0724 template<class T> std::vector<T> objects_sorted_by_values(const std::vector<T> & objects,
0725                           const std::vector<double> & values) {
0726   if (objects.size() != values.size()){
0727     throw Error("fjcore::objects_sorted_by_values(...): the size of the 'objects' vector must match the size of the 'values' vector");
0728   }
0729   std::vector<int> indices(values.size());
0730   for (size_t i = 0; i < indices.size(); i++) {indices[i] = i;}
0731   sort_indices(indices, values);
0732   std::vector<T> objects_sorted(objects.size());
0733   for (size_t i = 0; i < indices.size(); i++) {
0734     objects_sorted[i] = objects[indices[i]];
0735   }
0736   return objects_sorted;
0737 }
0738 class IndexedSortHelper {
0739 public:
0740   inline IndexedSortHelper (const std::vector<double> * reference_values) {
0741     _ref_values = reference_values;
0742   };
0743   inline int operator() (const int i1, const int i2) const {
0744     return  (*_ref_values)[i1] < (*_ref_values)[i2];
0745   };
0746 private:
0747   const std::vector<double> * _ref_values;
0748 };
0749 template <class L> inline  PseudoJet::PseudoJet(const L & some_four_vector) {
0750   reset(some_four_vector);
0751 }
0752 inline void PseudoJet::_reset_indices() {
0753   set_cluster_hist_index(-1);
0754   set_user_index(-1);
0755   _structure.reset();
0756   _user_info.reset();
0757 }
0758 inline double PseudoJet::m() const {
0759   double mm = m2();
0760   return mm < 0.0 ? -std::sqrt(-mm) : std::sqrt(mm);
0761 }
0762 inline void PseudoJet::reset(double px_in, double py_in, double pz_in, double E_in) {
0763   _px = px_in;
0764   _py = py_in;
0765   _pz = pz_in;
0766   _E  = E_in;
0767   _finish_init();
0768   _reset_indices();
0769 }
0770 inline void PseudoJet::reset_momentum(double px_in, double py_in, double pz_in, double E_in) {
0771   _px = px_in;
0772   _py = py_in;
0773   _pz = pz_in;
0774   _E  = E_in;
0775   _finish_init();
0776 }
0777 inline void PseudoJet::reset_momentum(const PseudoJet & pj) {
0778   _px  = pj._px ;
0779   _py  = pj._py ;
0780   _pz  = pj._pz ;
0781   _E   = pj._E  ;
0782   _phi = pj._phi;
0783   _rap = pj._rap;
0784   _kt2 = pj._kt2;
0785 }
0786 template<typename StructureType>
0787 const StructureType & PseudoJet::structure() const{
0788   return dynamic_cast<const StructureType &>(* validated_structure_ptr());
0789 }
0790 template<typename TransformerType>
0791 bool PseudoJet::has_structure_of() const{
0792   if (!_structure) return false;
0793   return dynamic_cast<const typename TransformerType::StructureType *>(_structure.get()) != 0;
0794 }
0795 template<typename TransformerType>
0796 const typename TransformerType::StructureType & PseudoJet::structure_of() const{
0797   if (!_structure)
0798     throw Error("Trying to access the structure of a PseudoJet without an associated structure");
0799   return dynamic_cast<const typename TransformerType::StructureType &>(*_structure);
0800 }
0801 PseudoJet join(const std::vector<PseudoJet> & pieces);
0802 PseudoJet join(const PseudoJet & j1);
0803 PseudoJet join(const PseudoJet & j1, const PseudoJet & j2);
0804 PseudoJet join(const PseudoJet & j1, const PseudoJet & j2, const PseudoJet & j3);
0805 PseudoJet join(const PseudoJet & j1, const PseudoJet & j2, const PseudoJet & j3, const PseudoJet & j4);
0806 FJCORE_END_NAMESPACE
0807 #endif // __FJCORE_PSEUDOJET_HH__
0808 #ifndef __FJCORE_FUNCTION_OF_PSEUDOJET_HH__
0809 #define __FJCORE_FUNCTION_OF_PSEUDOJET_HH__
0810 FJCORE_BEGIN_NAMESPACE
0811 template<typename TOut>
0812 class FunctionOfPseudoJet{
0813 public:
0814   FunctionOfPseudoJet(){}
0815   virtual ~FunctionOfPseudoJet(){}
0816   virtual std::string description() const{ return "";}
0817   virtual TOut result(const PseudoJet &pj) const = 0;
0818   TOut operator()(const PseudoJet &pj) const { return result(pj);}
0819   std::vector<TOut> operator()(const std::vector<PseudoJet> &pjs) const {
0820     std::vector<TOut> res(pjs.size());
0821     for (unsigned int i=0; i<pjs.size(); i++)
0822       res[i] = result(pjs[i]);
0823     return res;
0824   }
0825 };
0826 FJCORE_END_NAMESPACE
0827 #endif  // __FJCORE_FUNCTION_OF_PSEUDOJET_HH__
0828 #ifndef __FJCORE_SELECTOR_HH__
0829 #define __FJCORE_SELECTOR_HH__
0830 #include <limits>
0831 #include <cmath>
0832 FJCORE_BEGIN_NAMESPACE      // defined in fastjet/internal/base.hh
0833 class Selector;
0834 class SelectorWorker {
0835 public:
0836   virtual ~SelectorWorker() {}
0837   virtual bool pass(const PseudoJet & jet) const = 0;
0838   virtual void terminator(std::vector<const PseudoJet *> & jets) const {
0839     for (unsigned i = 0; i < jets.size(); i++) {
0840       if (jets[i] && !pass(*jets[i])) jets[i] = NULL;
0841     }
0842   }
0843   virtual bool applies_jet_by_jet() const {return true;}
0844   virtual std::string description() const {return "missing description";}
0845   virtual bool takes_reference() const { return false;}
0846   virtual void set_reference(const PseudoJet & /*reference*/){
0847     throw Error("set_reference(...) cannot be used for a selector worker that does not take a reference");
0848   }
0849   virtual SelectorWorker* copy(){
0850     throw Error("this SelectorWorker has nothing to copy");
0851   }
0852   virtual void get_rapidity_extent(double & rapmin, double & rapmax) const {
0853     rapmax = std::numeric_limits<double>::infinity();
0854     rapmin = -rapmax;
0855   }
0856   virtual bool is_geometric() const { return false;}
0857   virtual bool has_finite_area() const;
0858   virtual bool has_known_area() const { return false;}
0859   virtual double known_area() const{
0860     throw Error("this selector has no computable area");
0861   }
0862 };
0863 class Selector{
0864 public:
0865   Selector() {}
0866   Selector(SelectorWorker * worker_in) {_worker.reset(worker_in);}
0867   virtual ~Selector(){}
0868   bool pass(const PseudoJet & jet) const {
0869     if (!validated_worker()->applies_jet_by_jet()) {
0870       throw Error("Cannot apply this selector to an individual jet");
0871     }
0872     return _worker->pass(jet);
0873   }
0874   bool operator()(const PseudoJet & jet) const {
0875     return pass(jet);
0876   }
0877   unsigned int count(const std::vector<PseudoJet> & jets) const;
0878   PseudoJet sum(const std::vector<PseudoJet> & jets) const;
0879   double scalar_pt_sum(const std::vector<PseudoJet> & jets) const;
0880   void sift(const std::vector<PseudoJet> & jets,
0881           std::vector<PseudoJet> & jets_that_pass,
0882           std::vector<PseudoJet> & jets_that_fail) const;
0883   bool applies_jet_by_jet() const {
0884     return validated_worker()->applies_jet_by_jet();
0885   }
0886   std::vector<PseudoJet> operator()(const std::vector<PseudoJet> & jets) const;
0887   virtual void nullify_non_selected(std::vector<const PseudoJet *> & jets) const {
0888     validated_worker()->terminator(jets);
0889   }
0890   void get_rapidity_extent(double &rapmin, double &rapmax) const {
0891     return validated_worker()->get_rapidity_extent(rapmin, rapmax);
0892   }
0893   std::string description() const {
0894     return validated_worker()->description();
0895   }
0896   bool is_geometric() const{
0897     return validated_worker()->is_geometric();
0898   }
0899   bool has_finite_area() const{
0900     return validated_worker()->has_finite_area();
0901   }
0902   const SharedPtr<SelectorWorker> & worker() const {return _worker;}
0903   const SelectorWorker* validated_worker() const {
0904     const SelectorWorker* worker_ptr = _worker.get();
0905     if (worker_ptr == 0) throw InvalidWorker();
0906     return worker_ptr;
0907   }
0908   bool takes_reference() const {
0909     return validated_worker()->takes_reference();
0910   }
0911   const Selector & set_reference(const PseudoJet &reference){
0912     if (! validated_worker()->takes_reference()){
0913       return *this;
0914     }
0915     _copy_worker_if_needed();
0916     _worker->set_reference(reference);
0917     return *this;
0918   }
0919   class InvalidWorker : public Error {
0920   public:
0921     InvalidWorker() : Error("Attempt to use Selector with no valid underlying worker") {}
0922   };
0923   class InvalidArea : public Error {
0924   public:
0925     InvalidArea() : Error("Attempt to obtain area from Selector for which this is not meaningful") {}
0926   };
0927   Selector & operator &=(const Selector & b);
0928   Selector & operator |=(const Selector & b);
0929 protected:
0930   void _copy_worker_if_needed(){
0931     if (_worker.unique()) return;
0932     _worker.reset(_worker->copy());
0933   }
0934 private:
0935   SharedPtr<SelectorWorker> _worker; ///< the underlying worker
0936 };
0937 Selector SelectorIdentity();
0938 Selector operator!(const Selector & s);
0939 Selector operator ||(const Selector & s1, const Selector & s2);
0940 Selector operator&&(const Selector & s1, const Selector & s2);
0941 Selector operator*(const Selector & s1, const Selector & s2);
0942 Selector SelectorPtMin(double ptmin);                    ///< select objects with pt >= ptmin
0943 Selector SelectorPtMax(double ptmax);                    ///< select objects with pt <= ptmax
0944 Selector SelectorPtRange(double ptmin, double ptmax);    ///< select objects with ptmin <= pt <= ptmax
0945 Selector SelectorEtMin(double Etmin);                    ///< select objects with Et >= Etmin
0946 Selector SelectorEtMax(double Etmax);                    ///< select objects with Et <= Etmax
0947 Selector SelectorEtRange(double Etmin, double Etmax);    ///< select objects with Etmin <= Et <= Etmax
0948 Selector SelectorEMin(double Emin);                      ///< select objects with E >= Emin
0949 Selector SelectorEMax(double Emax);                      ///< select objects with E <= Emax
0950 Selector SelectorERange(double Emin, double Emax);       ///< select objects with Emin <= E <= Emax
0951 Selector SelectorMassMin(double Mmin);                      ///< select objects with Mass >= Mmin
0952 Selector SelectorMassMax(double Mmax);                      ///< select objects with Mass <= Mmax
0953 Selector SelectorMassRange(double Mmin, double Mmax);       ///< select objects with Mmin <= Mass <= Mmax
0954 Selector SelectorRapMin(double rapmin);                  ///< select objects with rap >= rapmin
0955 Selector SelectorRapMax(double rapmax);                  ///< select objects with rap <= rapmax
0956 Selector SelectorRapRange(double rapmin, double rapmax); ///< select objects with rapmin <= rap <= rapmax
0957 Selector SelectorAbsRapMin(double absrapmin);                     ///< select objects with |rap| >= absrapmin
0958 Selector SelectorAbsRapMax(double absrapmax);                     ///< select objects with |rap| <= absrapmax
0959 Selector SelectorAbsRapRange(double absrapmin, double absrapmax); ///< select objects with absrapmin <= |rap| <= absrapmax
0960 Selector SelectorEtaMin(double etamin);                  ///< select objects with eta >= etamin
0961 Selector SelectorEtaMax(double etamax);                  ///< select objects with eta <= etamax
0962 Selector SelectorEtaRange(double etamin, double etamax); ///< select objects with etamin <= eta <= etamax
0963 Selector SelectorAbsEtaMin(double absetamin);                     ///< select objects with |eta| >= absetamin
0964 Selector SelectorAbsEtaMax(double absetamax);                     ///< select objects with |eta| <= absetamax
0965 Selector SelectorAbsEtaRange(double absetamin, double absetamax); ///< select objects with absetamin <= |eta| <= absetamax
0966 Selector SelectorPhiRange(double phimin, double phimax); ///< select objects with phimin <= phi <= phimax
0967 Selector SelectorRapPhiRange(double rapmin, double rapmax, double phimin, double phimax);
0968 Selector SelectorNHardest(unsigned int n);
0969 Selector SelectorCircle(const double radius);
0970 Selector SelectorDoughnut(const double radius_in, const double radius_out);
0971 Selector SelectorStrip(const double half_width);
0972 Selector SelectorRectangle(const double half_rap_width, const double half_phi_width);
0973 Selector SelectorPtFractionMin(double fraction);
0974 Selector SelectorIsZero();
0975 FJCORE_END_NAMESPACE      // defined in fastjet/internal/base.hh
0976 #endif // __FJCORE_SELECTOR_HH__
0977 #ifndef __FJCORE_JETDEFINITION_HH__
0978 #define __FJCORE_JETDEFINITION_HH__
0979 #include<cassert>
0980 #include<string>
0981 #include<memory>
0982 FJCORE_BEGIN_NAMESPACE      // defined in fastjet/internal/base.hh
0983 std::string fastjet_version_string();
0984 enum Strategy {
0985   N2MHTLazy9AntiKtSeparateGhosts   = -10,
0986   N2MHTLazy9   = -7,
0987   N2MHTLazy25   = -6,
0988   N2MHTLazy9Alt   = -5,
0989   N2MinHeapTiled   = -4,
0990   N2Tiled     = -3,
0991   N2PoorTiled = -2,
0992   N2Plain     = -1,
0993   N3Dumb      =  0,
0994   Best        =  1,
0995   NlnN        =  2,
0996   NlnN3pi     =  3,
0997   NlnN4pi     =  4,
0998   NlnNCam4pi   = 14,
0999   NlnNCam2pi2R = 13,
1000   NlnNCam      = 12, // 2piMultD
1001   BestFJ30     =  21,
1002   plugin_strategy = 999
1003 };
1004 enum JetAlgorithm {
1005   kt_algorithm=0,
1006   cambridge_algorithm=1,
1007   antikt_algorithm=2,
1008   genkt_algorithm=3,
1009   cambridge_for_passive_algorithm=11,
1010   genkt_for_passive_algorithm=13,
1011   ee_kt_algorithm=50,
1012   ee_genkt_algorithm=53,
1013   plugin_algorithm = 99,
1014   undefined_jet_algorithm = 999
1015 };
1016 typedef JetAlgorithm JetFinder;
1017 const JetAlgorithm aachen_algorithm = cambridge_algorithm;
1018 const JetAlgorithm cambridge_aachen_algorithm = cambridge_algorithm;
1019 enum RecombinationScheme {
1020   E_scheme=0,
1021   pt_scheme=1,
1022   pt2_scheme=2,
1023   Et_scheme=3,
1024   Et2_scheme=4,
1025   BIpt_scheme=5,
1026   BIpt2_scheme=6,
1027   WTA_pt_scheme=7,
1028   WTA_modp_scheme=8,
1029   external_scheme = 99
1030 };
1031 class ClusterSequence;
1032 class JetDefinition {
1033 public:
1034   class Plugin;
1035   class Recombiner;
1036   JetDefinition(JetAlgorithm jet_algorithm_in,
1037                 double R_in,
1038                 RecombinationScheme recomb_scheme_in = E_scheme,
1039                 Strategy strategy_in = Best) {
1040     *this = JetDefinition(jet_algorithm_in, R_in, recomb_scheme_in, strategy_in, 1);
1041   }
1042   JetDefinition(JetAlgorithm jet_algorithm_in,
1043                 RecombinationScheme recomb_scheme_in = E_scheme,
1044                 Strategy strategy_in = Best) {
1045     double dummyR = 0.0;
1046     *this = JetDefinition(jet_algorithm_in, dummyR, recomb_scheme_in, strategy_in, 0);
1047   }
1048   JetDefinition(JetAlgorithm jet_algorithm_in,
1049                 double R_in,
1050                 double xtra_param_in,
1051                 RecombinationScheme recomb_scheme_in = E_scheme,
1052                 Strategy strategy_in = Best) {
1053     *this = JetDefinition(jet_algorithm_in, R_in, recomb_scheme_in, strategy_in, 2);
1054     set_extra_param(xtra_param_in);
1055   }
1056   JetDefinition(JetAlgorithm jet_algorithm_in,
1057                 double R_in,
1058                 const Recombiner * recombiner_in,
1059                 Strategy strategy_in = Best) {
1060     *this = JetDefinition(jet_algorithm_in, R_in, external_scheme, strategy_in);
1061     _recombiner = recombiner_in;
1062   }
1063   JetDefinition(JetAlgorithm jet_algorithm_in,
1064                 const Recombiner * recombiner_in,
1065                 Strategy strategy_in = Best) {
1066     *this = JetDefinition(jet_algorithm_in, external_scheme, strategy_in);
1067     _recombiner = recombiner_in;
1068   }
1069   JetDefinition(JetAlgorithm jet_algorithm_in,
1070                 double R_in,
1071                 double xtra_param_in,
1072                 const Recombiner * recombiner_in,
1073                 Strategy strategy_in = Best) {
1074     *this = JetDefinition(jet_algorithm_in, R_in, xtra_param_in, external_scheme, strategy_in);
1075     _recombiner = recombiner_in;
1076   }
1077   JetDefinition()  {
1078     *this = JetDefinition(undefined_jet_algorithm, 1.0);
1079   }
1080   JetDefinition(const Plugin * plugin_in) {
1081     _plugin = plugin_in;
1082     _strategy = plugin_strategy;
1083     _Rparam = _plugin->R();
1084     _extra_param = 0.0; // a dummy value to keep static code checkers happy
1085     _jet_algorithm = plugin_algorithm;
1086     set_recombination_scheme(E_scheme);
1087   }
1088   JetDefinition(JetAlgorithm jet_algorithm_in,
1089                 double R_in,
1090                 RecombinationScheme recomb_scheme_in,
1091                 Strategy strategy_in,
1092                 int nparameters_in);
1093   FJCORE_DEPRECATED_MSG("This argument ordering is deprecated. Use JetDefinition(alg, R, strategy, scheme[, n_parameters]) instead")
1094   JetDefinition(JetAlgorithm jet_algorithm_in,
1095                 double R_in,
1096                 Strategy strategy_in,
1097                 RecombinationScheme recomb_scheme_in = E_scheme,
1098                 int nparameters_in = 1){
1099     (*this) = JetDefinition(jet_algorithm_in,R_in,recomb_scheme_in,strategy_in,nparameters_in);
1100   }
1101   template <class L>
1102   std::vector<PseudoJet> operator()(const std::vector<L> & particles) const;
1103   static const double max_allowable_R; //= 1000.0;
1104   void set_recombination_scheme(RecombinationScheme);
1105   void set_recombiner(const Recombiner * recomb) {
1106     if (_shared_recombiner) _shared_recombiner.reset(recomb);
1107     _recombiner = recomb;
1108     _default_recombiner = DefaultRecombiner(external_scheme);
1109   }
1110   void set_recombiner(const JetDefinition &other_jet_def);
1111   void delete_recombiner_when_unused();
1112   const Plugin * plugin() const {return _plugin;};
1113   void delete_plugin_when_unused();
1114   JetAlgorithm jet_algorithm  () const {return _jet_algorithm  ;}
1115   JetAlgorithm jet_finder     () const {return _jet_algorithm  ;}
1116   double    R           () const {return _Rparam      ;}
1117   double    extra_param () const {return _extra_param ;}
1118   Strategy  strategy    () const {return _strategy    ;}
1119   RecombinationScheme recombination_scheme() const {
1120     return _default_recombiner.scheme();}
1121   void set_jet_algorithm(JetAlgorithm njf) {_jet_algorithm = njf;}
1122   void set_jet_finder(JetAlgorithm njf)    {_jet_algorithm = njf;}
1123   void set_extra_param(double xtra_param) {_extra_param = xtra_param;}
1124   const Recombiner * recombiner() const {
1125     return _recombiner == 0 ? & _default_recombiner : _recombiner;}
1126   bool has_same_recombiner(const JetDefinition &other_jd) const;
1127   bool is_spherical() const;
1128   std::string description() const;
1129   std::string description_no_recombiner() const;
1130   static std::string algorithm_description(const JetAlgorithm jet_alg);
1131   static unsigned int n_parameters_for_algorithm(const JetAlgorithm jet_alg);
1132 public:
1133   class Recombiner {
1134   public:
1135     virtual std::string description() const = 0;
1136     virtual void recombine(const PseudoJet & pa, const PseudoJet & pb,
1137                            PseudoJet & pab) const = 0;
1138     virtual void preprocess(PseudoJet & ) const {};
1139     virtual ~Recombiner() {};
1140     inline void plus_equal(PseudoJet & pa, const PseudoJet & pb) const {
1141       PseudoJet pres;
1142       recombine(pa,pb,pres);
1143       pa = pres;
1144     }
1145   };
1146   class DefaultRecombiner : public Recombiner {
1147   public:
1148     DefaultRecombiner(RecombinationScheme recomb_scheme = E_scheme) :
1149       _recomb_scheme(recomb_scheme) {}
1150     virtual std::string description() const FJCORE_OVERRIDE;
1151     virtual void recombine(const PseudoJet & pa, const PseudoJet & pb,
1152                            PseudoJet & pab) const FJCORE_OVERRIDE;
1153     virtual void preprocess(PseudoJet & p) const FJCORE_OVERRIDE;
1154     RecombinationScheme scheme() const {return _recomb_scheme;}
1155   private:
1156     RecombinationScheme _recomb_scheme;
1157   };
1158   class Plugin{
1159   public:
1160     virtual std::string description() const = 0;
1161     virtual void run_clustering(ClusterSequence &) const = 0;
1162     virtual double R() const = 0;
1163     virtual bool supports_ghosted_passive_areas() const {return false;}
1164     virtual void set_ghost_separation_scale(double scale) const;
1165     virtual double ghost_separation_scale() const {return 0.0;}
1166     virtual bool exclusive_sequence_meaningful() const {return false;}
1167     virtual bool is_spherical() const {return false;}
1168     virtual ~Plugin() {};
1169   };
1170 private:
1171   JetAlgorithm _jet_algorithm;
1172   double    _Rparam;
1173   double    _extra_param ; ///< parameter whose meaning varies according to context
1174   Strategy  _strategy  ;
1175   const Plugin * _plugin;
1176   SharedPtr<const Plugin> _plugin_shared;
1177   DefaultRecombiner _default_recombiner;
1178   const Recombiner * _recombiner;
1179   SharedPtr<const Recombiner> _shared_recombiner;
1180 };
1181 PseudoJet join(const std::vector<PseudoJet> & pieces, const JetDefinition::Recombiner & recombiner);
1182 PseudoJet join(const PseudoJet & j1,
1183            const JetDefinition::Recombiner & recombiner);
1184 PseudoJet join(const PseudoJet & j1, const PseudoJet & j2,
1185            const JetDefinition::Recombiner & recombiner);
1186 PseudoJet join(const PseudoJet & j1, const PseudoJet & j2, const PseudoJet & j3,
1187            const JetDefinition::Recombiner & recombiner);
1188 PseudoJet join(const PseudoJet & j1, const PseudoJet & j2, const PseudoJet & j3, const PseudoJet & j4,
1189            const JetDefinition::Recombiner & recombiner);
1190 FJCORE_END_NAMESPACE
1191 #endif // __FJCORE_JETDEFINITION_HH__
1192 #ifndef __FJCORE_COMPOSITEJET_STRUCTURE_HH__
1193 #define __FJCORE_COMPOSITEJET_STRUCTURE_HH__
1194 FJCORE_BEGIN_NAMESPACE      // defined in fastjet/internal/base.hh
1195 class CompositeJetStructure : public PseudoJetStructureBase{
1196 public:
1197   CompositeJetStructure() : _area_4vector_ptr(0){};
1198   CompositeJetStructure(const std::vector<PseudoJet> & initial_pieces,
1199             const JetDefinition::Recombiner * recombiner = 0);
1200   virtual ~CompositeJetStructure(){
1201     if (_area_4vector_ptr) delete _area_4vector_ptr;
1202   };
1203   virtual std::string description() const FJCORE_OVERRIDE;
1204   virtual bool has_constituents() const FJCORE_OVERRIDE;
1205   virtual std::vector<PseudoJet> constituents(const PseudoJet &jet) const FJCORE_OVERRIDE;
1206   virtual bool has_pieces(const PseudoJet & /*jet*/) const FJCORE_OVERRIDE {return true;}
1207   virtual std::vector<PseudoJet> pieces(const PseudoJet &jet) const FJCORE_OVERRIDE;
1208 protected:
1209   std::vector<PseudoJet> _pieces;  ///< the pieces building the jet
1210   PseudoJet * _area_4vector_ptr;   ///< pointer to the 4-vector jet area
1211 };
1212 template<typename T> PseudoJet join(const std::vector<PseudoJet> & pieces){
1213   PseudoJet result(0.0,0.0,0.0,0.0);
1214   for (unsigned int i=0; i<pieces.size(); i++){
1215     const PseudoJet it = pieces[i];
1216     result += it;
1217   }
1218   T *cj_struct = new T(pieces);
1219   result.set_structure_shared_ptr(SharedPtr<PseudoJetStructureBase>(cj_struct));
1220   return result;
1221 }
1222 template<typename T> PseudoJet join(const PseudoJet & j1){
1223   return join<T>(std::vector<PseudoJet>(1,j1));
1224 }
1225 template<typename T> PseudoJet join(const PseudoJet & j1, const PseudoJet & j2){
1226   std::vector<PseudoJet> pieces;
1227   pieces.push_back(j1);
1228   pieces.push_back(j2);
1229   return join<T>(pieces);
1230 }
1231 template<typename T> PseudoJet join(const PseudoJet & j1, const PseudoJet & j2,
1232                     const PseudoJet & j3){
1233   std::vector<PseudoJet> pieces;
1234   pieces.push_back(j1);
1235   pieces.push_back(j2);
1236   pieces.push_back(j3);
1237   return join<T>(pieces);
1238 }
1239 template<typename T> PseudoJet join(const PseudoJet & j1, const PseudoJet & j2,
1240                     const PseudoJet & j3, const PseudoJet & j4){
1241   std::vector<PseudoJet> pieces;
1242   pieces.push_back(j1);
1243   pieces.push_back(j2);
1244   pieces.push_back(j3);
1245   pieces.push_back(j4);
1246   return join<T>(pieces);
1247 }
1248 template<typename T> PseudoJet join(const std::vector<PseudoJet> & pieces,
1249                     const JetDefinition::Recombiner & recombiner){
1250   PseudoJet result;
1251   if (pieces.size()>0){
1252     result = pieces[0];
1253     for (unsigned int i=1; i<pieces.size(); i++){
1254       recombiner.plus_equal(result, pieces[i]);
1255     }
1256   }
1257   T *cj_struct = new T(pieces, &recombiner);
1258   result.set_structure_shared_ptr(SharedPtr<PseudoJetStructureBase>(cj_struct));
1259   return result;
1260 }
1261 template<typename T> PseudoJet join(const PseudoJet & j1,
1262                     const JetDefinition::Recombiner & recombiner){
1263   return join<T>(std::vector<PseudoJet>(1,j1), recombiner);
1264 }
1265 template<typename T> PseudoJet join(const PseudoJet & j1, const PseudoJet & j2,
1266                     const JetDefinition::Recombiner & recombiner){
1267   std::vector<PseudoJet> pieces;
1268   pieces.reserve(2);
1269   pieces.push_back(j1);
1270   pieces.push_back(j2);
1271   return join<T>(pieces, recombiner);
1272 }
1273 template<typename T> PseudoJet join(const PseudoJet & j1, const PseudoJet & j2,
1274                     const PseudoJet & j3,
1275                     const JetDefinition::Recombiner & recombiner){
1276   std::vector<PseudoJet> pieces;
1277   pieces.reserve(3);
1278   pieces.push_back(j1);
1279   pieces.push_back(j2);
1280   pieces.push_back(j3);
1281   return join<T>(pieces, recombiner);
1282 }
1283 template<typename T> PseudoJet join(const PseudoJet & j1, const PseudoJet & j2,
1284                     const PseudoJet & j3, const PseudoJet & j4,
1285                     const JetDefinition::Recombiner & recombiner){
1286   std::vector<PseudoJet> pieces;
1287   pieces.reserve(4);
1288   pieces.push_back(j1);
1289   pieces.push_back(j2);
1290   pieces.push_back(j3);
1291   pieces.push_back(j4);
1292   return join<T>(pieces, recombiner);
1293 }
1294 FJCORE_END_NAMESPACE      // defined in fastjet/internal/base.hh
1295 #endif // __FJCORE_MERGEDJET_STRUCTURE_HH__
1296 #ifndef __FJCORE_CLUSTER_SEQUENCE_STRUCTURE_HH__
1297 #define __FJCORE_CLUSTER_SEQUENCE_STRUCTURE_HH__
1298 #include <vector>
1299 FJCORE_BEGIN_NAMESPACE      // defined in fastjet/internal/base.hh
1300 class ClusterSequenceStructure : public PseudoJetStructureBase{
1301 public:
1302   ClusterSequenceStructure() : _associated_cs(NULL){}
1303   ClusterSequenceStructure(const ClusterSequence *cs){
1304     set_associated_cs(cs);
1305   };
1306   virtual ~ClusterSequenceStructure();
1307   virtual std::string description() const FJCORE_OVERRIDE{
1308     return "PseudoJet with an associated ClusterSequence";
1309   }
1310   virtual bool has_associated_cluster_sequence() const FJCORE_OVERRIDE{ return true;}
1311   virtual const ClusterSequence* associated_cluster_sequence() const FJCORE_OVERRIDE;
1312   virtual bool has_valid_cluster_sequence() const FJCORE_OVERRIDE;
1313   virtual const ClusterSequence * validated_cs() const FJCORE_OVERRIDE;
1314   virtual void set_associated_cs(const ClusterSequence * new_cs){
1315     _associated_cs = new_cs;
1316   }
1317   virtual bool has_partner(const PseudoJet &reference, PseudoJet &partner) const FJCORE_OVERRIDE;
1318   virtual bool has_child(const PseudoJet &reference, PseudoJet &child) const FJCORE_OVERRIDE;
1319   virtual bool has_parents(const PseudoJet &reference, PseudoJet &parent1, PseudoJet &parent2) const FJCORE_OVERRIDE;
1320   virtual bool object_in_jet(const PseudoJet &reference, const PseudoJet &jet) const FJCORE_OVERRIDE;
1321   virtual bool has_constituents() const FJCORE_OVERRIDE;
1322   virtual std::vector<PseudoJet> constituents(const PseudoJet &reference) const FJCORE_OVERRIDE;
1323   virtual bool has_exclusive_subjets() const FJCORE_OVERRIDE;
1324   virtual std::vector<PseudoJet> exclusive_subjets(const PseudoJet &reference, const double & dcut) const FJCORE_OVERRIDE;
1325   virtual int n_exclusive_subjets(const PseudoJet &reference, const double & dcut) const FJCORE_OVERRIDE;
1326   virtual std::vector<PseudoJet> exclusive_subjets_up_to (const PseudoJet &reference, int nsub) const FJCORE_OVERRIDE;
1327   virtual double exclusive_subdmerge(const PseudoJet &reference, int nsub) const FJCORE_OVERRIDE;
1328   virtual double exclusive_subdmerge_max(const PseudoJet &reference, int nsub) const FJCORE_OVERRIDE;
1329   virtual bool has_pieces(const PseudoJet &reference) const FJCORE_OVERRIDE;
1330   virtual std::vector<PseudoJet> pieces(const PseudoJet &reference) const FJCORE_OVERRIDE;
1331 protected:
1332   const ClusterSequence *_associated_cs;
1333 };
1334 FJCORE_END_NAMESPACE
1335 #endif  //  __FJCORE_CLUSTER_SEQUENCE_STRUCTURE_HH__
1336 #ifndef __FJCORE_CLUSTERSEQUENCE_HH__
1337 #define __FJCORE_CLUSTERSEQUENCE_HH__
1338 #include<vector>
1339 #include<map>
1340 #include<memory>
1341 #include<cassert>
1342 #include<iostream>
1343 #include<string>
1344 #include<set>
1345 #include<cmath> // needed to get double std::abs(double)
1346 FJCORE_BEGIN_NAMESPACE      // defined in fastjet/internal/base.hh
1347 class ClusterSequenceStructure;
1348 class DynamicNearestNeighbours;
1349 class ClusterSequence {
1350  public:
1351   ClusterSequence () : _deletes_self_when_unused(false) {}
1352   template<class L> ClusterSequence (
1353                       const std::vector<L> & pseudojets,
1354                   const JetDefinition & jet_def,
1355                   const bool & writeout_combinations = false);
1356   ClusterSequence (const ClusterSequence & cs) : _deletes_self_when_unused(false) {
1357     transfer_from_sequence(cs);
1358   }
1359   ClusterSequence & operator=(const ClusterSequence & cs);
1360   virtual ~ClusterSequence (); //{}
1361   std::vector<PseudoJet> inclusive_jets (const double ptmin = 0.0) const;
1362   int n_exclusive_jets (const double dcut) const;
1363   std::vector<PseudoJet> exclusive_jets (const double dcut) const;
1364   std::vector<PseudoJet> exclusive_jets (const int njets) const;
1365   std::vector<PseudoJet> exclusive_jets_up_to (const int njets) const;
1366   double exclusive_dmerge (const int njets) const;
1367   double exclusive_dmerge_max (const int njets) const;
1368   double exclusive_ymerge (int njets) const {return exclusive_dmerge(njets) / Q2();}
1369   double exclusive_ymerge_max (int njets) const {return exclusive_dmerge_max(njets)/Q2();}
1370   int n_exclusive_jets_ycut (double ycut) const {return n_exclusive_jets(ycut*Q2());}
1371   std::vector<PseudoJet> exclusive_jets_ycut (double ycut) const {
1372     int njets = n_exclusive_jets_ycut(ycut);
1373     return exclusive_jets(njets);
1374   }
1375   std::vector<PseudoJet> exclusive_subjets (const PseudoJet & jet,
1376                                             const double dcut) const;
1377   int n_exclusive_subjets(const PseudoJet & jet,
1378                           const double dcut) const;
1379   std::vector<PseudoJet> exclusive_subjets (const PseudoJet & jet,
1380                                             int nsub) const;
1381   std::vector<PseudoJet> exclusive_subjets_up_to (const PseudoJet & jet,
1382                           int nsub) const;
1383   double exclusive_subdmerge(const PseudoJet & jet, int nsub) const;
1384   double exclusive_subdmerge_max(const PseudoJet & jet, int nsub) const;
1385   double Q() const {return _Qtot;}
1386   double Q2() const {return _Qtot*_Qtot;}
1387   bool object_in_jet(const PseudoJet & object, const PseudoJet & jet) const;
1388   bool has_parents(const PseudoJet & jet, PseudoJet & parent1,
1389                PseudoJet & parent2) const;
1390   bool has_child(const PseudoJet & jet, PseudoJet & child) const;
1391   bool has_child(const PseudoJet & jet, const PseudoJet * & childp) const;
1392   bool has_partner(const PseudoJet & jet, PseudoJet & partner) const;
1393   std::vector<PseudoJet> constituents (const PseudoJet & jet) const;
1394   void print_jets_for_root(const std::vector<PseudoJet> & jets,
1395                            std::ostream & ostr = std::cout) const;
1396   void print_jets_for_root(const std::vector<PseudoJet> & jets,
1397                            const std::string & filename,
1398                const std::string & comment = "") const;
1399   void add_constituents (const PseudoJet & jet,
1400              std::vector<PseudoJet> & subjet_vector) const;
1401   inline Strategy strategy_used () const {return _strategy;}
1402   std::string strategy_string () const {return strategy_string(_strategy);}
1403   std::string strategy_string (Strategy strategy_in) const;
1404   const JetDefinition & jet_def() const {return _jet_def;}
1405   void delete_self_when_unused();
1406   bool will_delete_self_when_unused() const {return _deletes_self_when_unused;}
1407   void signal_imminent_self_deletion() const;
1408   double jet_scale_for_algorithm(const PseudoJet & jet) const;
1409   void plugin_record_ij_recombination(int jet_i, int jet_j, double dij,
1410                       int & newjet_k) {
1411     assert(plugin_activated());
1412     _do_ij_recombination_step(jet_i, jet_j, dij, newjet_k);
1413   }
1414   void plugin_record_ij_recombination(int jet_i, int jet_j, double dij,
1415                       const PseudoJet & newjet,
1416                       int & newjet_k);
1417   void plugin_record_iB_recombination(int jet_i, double diB) {
1418     assert(plugin_activated());
1419     _do_iB_recombination_step(jet_i, diB);
1420   }
1421   class Extras {
1422   public:
1423     virtual ~Extras() {}
1424     virtual std::string description() const {return "This is a dummy extras class that contains no extra information! Derive from it if you want to use it to provide extra information from a plugin jet finder";}
1425   };
1426   inline void plugin_associate_extras(Extras * extras_in) {
1427     _extras.reset(extras_in);
1428   }
1429 #ifdef FJCORE_HAVE_AUTO_PTR_INTERFACE
1430   FJCORE_DEPRECATED_MSG("Please use ClusterSequence::plugin_associate_extras(Extras * extras_in)) instead")
1431   inline void plugin_associate_extras(std::auto_ptr<Extras> extras_in){
1432     _extras.reset(extras_in.release());
1433   }
1434 #endif
1435   inline bool plugin_activated() const {return _plugin_activated;}
1436   const Extras * extras() const {return _extras.get();}
1437   template<class GBJ> void plugin_simple_N2_cluster () {
1438     assert(plugin_activated());
1439     _simple_N2_cluster<GBJ>();
1440   }
1441 public:
1442   struct history_element{
1443     int parent1; /// index in _history where first parent of this jet
1444     int parent2; /// index in _history where second parent of this jet
1445     int child;   /// index in _history where the current jet is
1446          /// recombined with another jet to form its child. It
1447          /// is Invalid if this jet does not further
1448          /// recombine.
1449     int jetp_index; /// index in the _jets vector where we will find the
1450     double dij;  /// the distance corresponding to the recombination
1451          /// at this stage of the clustering.
1452     double max_dij_so_far; /// the largest recombination distance seen
1453                /// so far in the clustering history.
1454   };
1455   enum JetType {Invalid=-3, InexistentParent = -2, BeamJet = -1};
1456   const std::vector<PseudoJet> & jets()    const;
1457   const std::vector<history_element> & history() const;
1458   unsigned int n_particles() const;
1459   std::vector<int> particle_jet_indices(const std::vector<PseudoJet> &) const;
1460   std::vector<int> unique_history_order() const;
1461   std::vector<PseudoJet> unclustered_particles() const;
1462   std::vector<PseudoJet> childless_pseudojets() const;
1463   bool contains(const PseudoJet & object) const;
1464   void transfer_from_sequence(const ClusterSequence & from_seq,
1465                   const FunctionOfPseudoJet<PseudoJet> * action_on_jets = 0);
1466   const SharedPtr<PseudoJetStructureBase> & structure_shared_ptr() const{
1467     return _structure_shared_ptr;
1468   }
1469   typedef ClusterSequenceStructure StructureType;
1470   static void print_banner();
1471   static void set_fastjet_banner_stream(std::ostream * ostr) {_fastjet_banner_ostr = ostr;}
1472   static std::ostream * fastjet_banner_stream() {return _fastjet_banner_ostr;}
1473 private:
1474   static std::ostream * _fastjet_banner_ostr;
1475 protected:
1476   JetDefinition _jet_def;
1477   template<class L> void _transfer_input_jets(
1478                                      const std::vector<L> & pseudojets);
1479   void _initialise_and_run (const JetDefinition & jet_def,
1480                 const bool & writeout_combinations);
1481   void _initialise_and_run_no_decant();
1482   void _decant_options(const JetDefinition & jet_def,
1483                        const bool & writeout_combinations);
1484   void _decant_options_partial();
1485   void _fill_initial_history();
1486   void _do_ij_recombination_step(const int jet_i, const int jet_j,
1487                  const double dij, int & newjet_k);
1488   void _do_iB_recombination_step(const int jet_i, const double diB);
1489   void _set_structure_shared_ptr(PseudoJet & j);
1490   void _update_structure_use_count();
1491   Strategy _best_strategy() const;
1492   class _Parabola {
1493   public:
1494     _Parabola(double a, double b, double c) : _a(a), _b(b), _c(c) {}
1495     inline double operator()(const double R) const {return _c*(_a*R*R + _b*R + 1);}
1496   private:
1497     double _a, _b, _c;
1498   };
1499   class _Line {
1500   public:
1501     _Line(double a, double b) : _a(a), _b(b) {}
1502     inline double operator()(const double R) const {return _a*R + _b;}
1503   private:
1504     double _a, _b;
1505   };
1506   std::vector<PseudoJet> _jets;
1507   std::vector<history_element> _history;
1508   void get_subhist_set(std::set<const history_element*> & subhist,
1509                        const  PseudoJet & jet, double dcut, int maxjet) const;
1510   bool _writeout_combinations;
1511   int  _initial_n;
1512   double _Rparam, _R2, _invR2;
1513   double _Qtot;
1514   Strategy    _strategy;
1515   JetAlgorithm  _jet_algorithm;
1516   SharedPtr<PseudoJetStructureBase> _structure_shared_ptr; //< will actually be of type ClusterSequenceStructure
1517   int _structure_use_count_after_construction; //< info of use when CS handles its own memory
1518   mutable bool _deletes_self_when_unused;
1519  private:
1520   bool _plugin_activated;
1521   SharedPtr<Extras> _extras; // things the plugin might want to add
1522   void _really_dumb_cluster ();
1523   void _delaunay_cluster ();
1524   template<class BJ> void _simple_N2_cluster ();
1525   void _tiled_N2_cluster ();
1526   void _faster_tiled_N2_cluster ();
1527   void _minheap_faster_tiled_N2_cluster();
1528   void _CP2DChan_cluster();
1529   void _CP2DChan_cluster_2pi2R ();
1530   void _CP2DChan_cluster_2piMultD ();
1531   void _CP2DChan_limited_cluster(double D);
1532   void _do_Cambridge_inclusive_jets();
1533   void _fast_NsqrtN_cluster();
1534   void _add_step_to_history( //const int step_number,
1535                             const int parent1,
1536                 const int parent2, const int jetp_index,
1537                 const double dij);
1538   void _extract_tree_children(int pos, std::valarray<bool> &,
1539         const std::valarray<int> &, std::vector<int> &) const;
1540   void _extract_tree_parents (int pos, std::valarray<bool> &,
1541                 const std::valarray<int> &,  std::vector<int> &) const;
1542   typedef std::pair<int,int> TwoVertices;
1543   typedef std::pair<double,TwoVertices> DijEntry;
1544   typedef std::multimap<double,TwoVertices> DistMap;
1545   void _add_ktdistance_to_map(const int ii,
1546                   DistMap & DijMap,
1547                   const DynamicNearestNeighbours * DNN);
1548   static bool _first_time;
1549   static LimitedWarning _exclusive_warnings;
1550   static LimitedWarning _changed_strategy_warning;
1551   struct BriefJet {
1552     double     eta, phi, kt2, NN_dist;
1553     BriefJet * NN;
1554     int        _jets_index;
1555   };
1556   class TiledJet {
1557   public:
1558     double     eta, phi, kt2, NN_dist;
1559     TiledJet * NN, *previous, * next;
1560     int        _jets_index, tile_index, diJ_posn;
1561     inline void label_minheap_update_needed() {diJ_posn = 1;}
1562     inline void label_minheap_update_done()   {diJ_posn = 0;}
1563     inline bool minheap_update_needed() const {return diJ_posn==1;}
1564   };
1565   template <class J> void _bj_set_jetinfo( J * const jet,
1566                          const int _jets_index) const;
1567   void _bj_remove_from_tiles( TiledJet * const jet) const;
1568   template <class J> double _bj_dist(const J * const jeta,
1569             const J * const jetb) const;
1570   template <class J> double _bj_diJ(const J * const jeta) const;
1571   template <class J> inline J * _bj_of_hindex(
1572                           const int hist_index,
1573               J * const head, J * const tail)
1574     const {
1575     J * res;
1576     for(res = head; res<tail; res++) {
1577       if (_jets[res->_jets_index].cluster_hist_index() == hist_index) {break;}
1578     }
1579     return res;
1580   }
1581   template <class J> void _bj_set_NN_nocross(J * const jeta,
1582             J * const head, const J * const tail) const;
1583   template <class J> void _bj_set_NN_crosscheck(J * const jeta,
1584             J * const head, const J * const tail) const;
1585   static const int n_tile_neighbours = 9;
1586   struct Tile {
1587     Tile *   begin_tiles[n_tile_neighbours];
1588     Tile **  surrounding_tiles;
1589     Tile **  RH_tiles;
1590     Tile **  end_tiles;
1591     TiledJet * head;
1592     bool     tagged;
1593   };
1594   std::vector<Tile> _tiles;
1595   double _tiles_eta_min, _tiles_eta_max;
1596   double _tile_size_eta, _tile_size_phi;
1597   int    _n_tiles_phi,_tiles_ieta_min,_tiles_ieta_max;
1598   inline int _tile_index (int ieta, int iphi) const {
1599     return (ieta-_tiles_ieta_min)*_n_tiles_phi
1600                   + (iphi+_n_tiles_phi) % _n_tiles_phi;
1601   }
1602   int  _tile_index(const double eta, const double phi) const;
1603   void _tj_set_jetinfo ( TiledJet * const jet, const int _jets_index);
1604   void  _bj_remove_from_tiles(TiledJet * const jet);
1605   void _initialise_tiles();
1606   void _print_tiles(TiledJet * briefjets ) const;
1607   void _add_neighbours_to_tile_union(const int tile_index,
1608          std::vector<int> & tile_union, int & n_near_tiles) const;
1609   void _add_untagged_neighbours_to_tile_union(const int tile_index,
1610          std::vector<int> & tile_union, int & n_near_tiles);
1611   struct EEBriefJet {
1612     double NN_dist;  // obligatorily present
1613     double kt2;      // obligatorily present == E^2 in general
1614     EEBriefJet * NN; // must be present too
1615     int    _jets_index; // must also be present!
1616     double nx, ny, nz;  // our internal storage for fast distance calcs
1617   };
1618   void _simple_N2_cluster_BriefJet();
1619   void _simple_N2_cluster_EEBriefJet();
1620 };
1621 template<class L> void ClusterSequence::_transfer_input_jets(
1622                                        const std::vector<L> & pseudojets) {
1623   _jets.reserve(pseudojets.size()*2);
1624   for (unsigned int i = 0; i < pseudojets.size(); i++) {
1625     _jets.push_back(pseudojets[i]);}
1626 }
1627 template<class L> ClusterSequence::ClusterSequence (
1628                       const std::vector<L> & pseudojets,
1629                   const JetDefinition & jet_def_in,
1630                   const bool & writeout_combinations) :
1631   _jet_def(jet_def_in), _writeout_combinations(writeout_combinations),
1632   _structure_shared_ptr(new ClusterSequenceStructure(this))
1633 {
1634   _transfer_input_jets(pseudojets);
1635   _decant_options_partial();
1636   _initialise_and_run_no_decant();
1637 }
1638 inline const std::vector<PseudoJet> & ClusterSequence::jets () const {
1639   return _jets;
1640 }
1641 inline const std::vector<ClusterSequence::history_element> & ClusterSequence::history () const {
1642   return _history;
1643 }
1644 inline unsigned int ClusterSequence::n_particles() const {return _initial_n;}
1645 #ifndef __CINT__
1646 template<class L>
1647 std::vector<PseudoJet> JetDefinition::operator()(const std::vector<L> & particles) const {
1648   ClusterSequence * cs = new ClusterSequence(particles, *this);
1649   std::vector<PseudoJet> jets;
1650   if (is_spherical()) {
1651     jets = sorted_by_E(cs->inclusive_jets());
1652   } else {
1653     jets = sorted_by_pt(cs->inclusive_jets());
1654   }
1655   if (jets.size() != 0) {
1656     cs->delete_self_when_unused();
1657   } else {
1658     delete cs;
1659   }
1660   return jets;
1661 }
1662 #endif // __CINT__
1663 template <class J> inline void ClusterSequence::_bj_set_jetinfo(
1664                             J * const jetA, const int _jets_index) const {
1665     jetA->eta  = _jets[_jets_index].rap();
1666     jetA->phi  = _jets[_jets_index].phi_02pi();
1667     jetA->kt2  = jet_scale_for_algorithm(_jets[_jets_index]);
1668     jetA->_jets_index = _jets_index;
1669     jetA->NN_dist = _R2;
1670     jetA->NN      = NULL;
1671 }
1672 template <class J> inline double ClusterSequence::_bj_dist(
1673                 const J * const jetA, const J * const jetB) const {
1674 #ifndef FJCORE_NEW_DELTA_PHI
1675   double dphi = std::abs(jetA->phi - jetB->phi);
1676   double deta = (jetA->eta - jetB->eta);
1677   if (dphi > pi) {dphi = twopi - dphi;}
1678 #else
1679   double dphi = pi-std::abs(pi-std::abs(jetA->phi - jetB->phi));
1680   double deta = (jetA->eta - jetB->eta);
1681 #endif
1682   return dphi*dphi + deta*deta;
1683 }
1684 template <class J> inline double ClusterSequence::_bj_diJ(const J * const jet) const {
1685   double kt2 = jet->kt2;
1686   if (jet->NN != NULL) {if (jet->NN->kt2 < kt2) {kt2 = jet->NN->kt2;}}
1687   return jet->NN_dist * kt2;
1688 }
1689 template <class J> inline void ClusterSequence::_bj_set_NN_nocross(
1690                  J * const jet, J * const head, const J * const tail) const {
1691   double NN_dist = _R2;
1692   J * NN  = NULL;
1693   if (head < jet) {
1694     for (J * jetB = head; jetB != jet; jetB++) {
1695       double dist = _bj_dist(jet,jetB);
1696       if (dist < NN_dist) {
1697     NN_dist = dist;
1698     NN = jetB;
1699       }
1700     }
1701   }
1702   if (tail > jet) {
1703     for (J * jetB = jet+1; jetB != tail; jetB++) {
1704       double dist = _bj_dist(jet,jetB);
1705       if (dist < NN_dist) {
1706     NN_dist = dist;
1707     NN = jetB;
1708       }
1709     }
1710   }
1711   jet->NN = NN;
1712   jet->NN_dist = NN_dist;
1713 }
1714 template <class J> inline void ClusterSequence::_bj_set_NN_crosscheck(J * const jet,
1715             J * const head, const J * const tail) const {
1716   double NN_dist = _R2;
1717   J * NN  = NULL;
1718   for (J * jetB = head; jetB != tail; jetB++) {
1719     double dist = _bj_dist(jet,jetB);
1720     if (dist < NN_dist) {
1721       NN_dist = dist;
1722       NN = jetB;
1723     }
1724     if (dist < jetB->NN_dist) {
1725       jetB->NN_dist = dist;
1726       jetB->NN = jet;
1727     }
1728   }
1729   jet->NN = NN;
1730   jet->NN_dist = NN_dist;
1731 }
1732 FJCORE_END_NAMESPACE
1733 #endif // __FJCORE_CLUSTERSEQUENCE_HH__
1734 #ifndef __FJCORE_NNBASE_HH__
1735 #define __FJCORE_NNBASE_HH__
1736 FJCORE_BEGIN_NAMESPACE      // defined in fastjet/internal/base.hh
1737 class _NoInfo {};
1738 template<class I> class NNInfo {
1739 public:
1740   NNInfo()         : _info(NULL) {}
1741   NNInfo(I * info) : _info(info) {}
1742   template<class BJ> void init_jet(BJ * briefjet, const fjcore::PseudoJet & jet, int index) { briefjet->init(jet, index, _info);}
1743 private:
1744   I * _info;
1745 };
1746 template<> class NNInfo<_NoInfo>  {
1747 public:
1748   NNInfo()           {}
1749   NNInfo(_NoInfo * ) {}
1750   template<class BJ> void init_jet(BJ * briefjet, const fjcore::PseudoJet & jet, int index) { briefjet->init(jet, index);}
1751 };
1752 template<class I = _NoInfo> class NNBase : public NNInfo<I> {
1753 public:
1754   NNBase() {}
1755   NNBase(I * info) : NNInfo<I>(info) {}
1756   virtual void start(const std::vector<PseudoJet> & jets) = 0;
1757   virtual double dij_min(int & iA, int & iB) = 0;
1758   virtual void remove_jet(int iA) = 0;
1759   virtual void merge_jets(int iA, int iB, const PseudoJet & jet, int jet_index) =  0;
1760   virtual ~NNBase() {};
1761 };
1762 FJCORE_END_NAMESPACE      // defined in fastjet/internal/base.hh
1763 #endif // __FJCORE_NNBASE_HH__
1764 #ifndef __FJCORE_NNH_HH__
1765 #define __FJCORE_NNH_HH__
1766 FJCORE_BEGIN_NAMESPACE      // defined in fastjet/internal/base.hh
1767 template<class BJ, class I = _NoInfo> class NNH : public NNBase<I> {
1768 public:
1769   NNH(const std::vector<PseudoJet> & jets)           : NNBase<I>()     {start(jets);}
1770   NNH(const std::vector<PseudoJet> & jets, I * info) : NNBase<I>(info) {start(jets);}
1771   void start(const std::vector<PseudoJet> & jets);
1772   double dij_min(int & iA, int & iB);
1773   void remove_jet(int iA);
1774   void merge_jets(int iA, int iB, const PseudoJet & jet, int jet_index);
1775   ~NNH() {
1776     delete[] briefjets;
1777   }
1778 private:
1779   class NNBJ; // forward declaration
1780   void set_NN_crosscheck(NNBJ * jet, NNBJ * begin, NNBJ * end);
1781   void set_NN_nocross   (NNBJ * jet, NNBJ * begin, NNBJ * end);
1782   NNBJ * briefjets;
1783   NNBJ * head, * tail;
1784   int n;
1785   std::vector<NNBJ *> where_is;
1786   class NNBJ : public BJ {
1787   public:
1788     void init(const PseudoJet & jet, int index_in) {
1789       BJ::init(jet);
1790       other_init(index_in);
1791     }
1792     void init(const PseudoJet & jet, int index_in, I * info) {
1793       BJ::init(jet, info);
1794       other_init(index_in);
1795     }
1796     void other_init(int index_in) {
1797       _index = index_in;
1798       NN_dist = BJ::beam_distance();
1799       NN = NULL;
1800     }
1801     int index() const {return _index;}
1802     double NN_dist;
1803     NNBJ * NN;
1804   private:
1805     int _index;
1806   };
1807 };
1808 template<class BJ, class I> void NNH<BJ,I>::start(const std::vector<PseudoJet> & jets) {
1809   n = jets.size();
1810   briefjets = new NNBJ[n];
1811   where_is.resize(2*n);
1812   NNBJ * jetA = briefjets;
1813   for (int i = 0; i< n; i++) {
1814     this->init_jet(jetA, jets[i], i);
1815     where_is[i] = jetA;
1816     jetA++; // move on to next entry of briefjets
1817   }
1818   tail = jetA; // a semaphore for the end of briefjets
1819   head = briefjets; // a nicer way of naming start
1820   for (jetA = head + 1; jetA != tail; jetA++) {
1821     set_NN_crosscheck(jetA, head, jetA);
1822   }
1823 }
1824 template<class BJ, class I> double NNH<BJ,I>::dij_min(int & iA, int & iB) {
1825   double diJ_min = briefjets[0].NN_dist;
1826   int diJ_min_jet = 0;
1827   for (int i = 1; i < n; i++) {
1828     if (briefjets[i].NN_dist < diJ_min) {
1829       diJ_min_jet = i;
1830       diJ_min  = briefjets[i].NN_dist;
1831     }
1832   }
1833   NNBJ * jetA = & briefjets[diJ_min_jet];
1834   iA = jetA->index();
1835   iB = jetA->NN ? jetA->NN->index() : -1;
1836   return diJ_min;
1837 }
1838 template<class BJ, class I> void NNH<BJ,I>::remove_jet(int iA) {
1839   NNBJ * jetA = where_is[iA];
1840   tail--; n--;
1841   *jetA = *tail;
1842   where_is[jetA->index()] = jetA;
1843   for (NNBJ * jetI = head; jetI != tail; jetI++) {
1844     if (jetI->NN == jetA) set_NN_nocross(jetI, head, tail);
1845     if (jetI->NN == tail) {jetI->NN = jetA;}
1846   }
1847 }
1848 template<class BJ, class I> void NNH<BJ,I>::merge_jets(int iA, int iB,
1849                     const PseudoJet & jet, int index) {
1850   NNBJ * jetA = where_is[iA];
1851   NNBJ * jetB = where_is[iB];
1852   if (jetA < jetB) std::swap(jetA,jetB);
1853   this->init_jet(jetB, jet, index);
1854   if (index >= int(where_is.size())) where_is.resize(2*index);
1855   where_is[jetB->index()] = jetB;
1856   tail--; n--;
1857   *jetA = *tail;
1858   where_is[jetA->index()] = jetA;
1859   for (NNBJ * jetI = head; jetI != tail; jetI++) {
1860     if (jetI->NN == jetA || jetI->NN == jetB) {
1861       set_NN_nocross(jetI, head, tail);
1862     }
1863     double dist = jetI->distance(jetB);
1864     if (dist < jetI->NN_dist) {
1865       if (jetI != jetB) {
1866     jetI->NN_dist = dist;
1867     jetI->NN = jetB;
1868       }
1869     }
1870     if (dist < jetB->NN_dist) {
1871       if (jetI != jetB) {
1872     jetB->NN_dist = dist;
1873     jetB->NN      = jetI;
1874       }
1875     }
1876     if (jetI->NN == tail) {jetI->NN = jetA;}
1877   }
1878 }
1879 template <class BJ, class I> void NNH<BJ,I>::set_NN_crosscheck(NNBJ * jet,
1880             NNBJ * begin, NNBJ * end) {
1881   double NN_dist = jet->beam_distance();
1882   NNBJ * NN      = NULL;
1883   for (NNBJ * jetB = begin; jetB != end; jetB++) {
1884     double dist = jet->distance(jetB);
1885     if (dist < NN_dist) {
1886       NN_dist = dist;
1887       NN = jetB;
1888     }
1889     if (dist < jetB->NN_dist) {
1890       jetB->NN_dist = dist;
1891       jetB->NN = jet;
1892     }
1893   }
1894   jet->NN = NN;
1895   jet->NN_dist = NN_dist;
1896 }
1897 template <class BJ, class I>  void NNH<BJ,I>::set_NN_nocross(
1898                  NNBJ * jet, NNBJ * begin, NNBJ * end) {
1899   double NN_dist = jet->beam_distance();
1900   NNBJ * NN      = NULL;
1901   if (begin < jet) {
1902     for (NNBJ * jetB = begin; jetB != jet; jetB++) {
1903       double dist = jet->distance(jetB);
1904       if (dist < NN_dist) {
1905     NN_dist = dist;
1906     NN = jetB;
1907       }
1908     }
1909   }
1910   if (end > jet) {
1911     for (NNBJ * jetB = jet+1; jetB != end; jetB++) {
1912       double dist = jet->distance (jetB);
1913       if (dist < NN_dist) {
1914     NN_dist = dist;
1915     NN = jetB;
1916       }
1917     }
1918   }
1919   jet->NN = NN;
1920   jet->NN_dist = NN_dist;
1921 }
1922 FJCORE_END_NAMESPACE      // defined in fastjet/internal/base.hh
1923 #endif // __FJCORE_NNH_HH__
1924 #endif