Back to home page

EIC code displayed by LXR
 
 

    

File indexing completed on 2025-04-19 09:09:56 

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

This page was automatically generated by the 2.3.7 LXR engine.

The LXR team
Valid CSS 2.1! Valid HTML 4.01!