StRoot  1
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Groups Pages
FJcore.h
1 // fjcore -- extracted from FastJet v3.2.1 (http://fastjet.fr)
2 //
3 // fjcore constitutes a digest of the main FastJet functionality.
4 // The files fjcore.hh and fjcore.cc are meant to provide easy access to these
5 // core functions, in the form of single files and without the need of a full
6 // FastJet installation:
7 //
8 // g++ main.cc fjcore.cc
9 //
10 // with main.cc including fjcore.hh.
11 //
12 // A fortran interface, fjcorefortran.cc, is also provided. See the example
13 // and the Makefile for instructions.
14 //
15 // The results are expected to be identical to those obtained by linking to
16 // the full FastJet distribution.
17 //
18 // NOTE THAT, IN ORDER TO MAKE IT POSSIBLE FOR FJCORE AND THE FULL FASTJET
19 // TO COEXIST, THE FORMER USES THE "fjcore" NAMESPACE INSTEAD OF "fastjet".
20 //
21 // In particular, fjcore provides:
22 //
23 // - access to all native pp and ee algorithms, kt, anti-kt, C/A.
24 // For C/A, the NlnN method is available, while anti-kt and kt
25 // are limited to the N^2 one (still the fastest for N < 100k particles)
26 // - access to selectors, for implementing cuts and selections
27 // - access to all functionalities related to pseudojets (e.g. a jet's
28 // structure or user-defined information)
29 //
30 // Instead, it does NOT provide:
31 //
32 // - jet areas functionality
33 // - background estimation
34 // - access to other algorithms via plugins
35 // - interface to CGAL
36 // - fastjet tools, e.g. filters, taggers
37 //
38 // If these functionalities are needed, the full FastJet installation must be
39 // used. The code will be fully compatible, with the sole replacement of the
40 // header files and of the fjcore namespace with the fastjet one.
41 //
42 // fjcore.hh and fjcore.cc are not meant to be human-readable.
43 // For documentation, see the full FastJet manual and doxygen at http://fastjet.fr
44 //
45 // Like FastJet, fjcore is released under the terms of the GNU General Public
46 // License version 2 (GPLv2). If you use this code as part of work towards a
47 // scientific publication, whether directly or contained within another program
48 // (e.g. Delphes, MadGraph, SpartyJet, Rivet, LHC collaboration software frameworks,
49 // etc.), you should include a citation to
50 //
51 // EPJC72(2012)1896 [arXiv:1111.6097] (FastJet User Manual)
52 // and, optionally, Phys.Lett.B641 (2006) 57 [arXiv:hep-ph/0512210]
53 //
54 // Copyright (c) 2005-2016, Matteo Cacciari, Gavin P. Salam and Gregory Soyez
55 //
56 //----------------------------------------------------------------------
57 // This file is part of FastJet.
58 //
59 // FastJet is free software; you can redistribute it and/or modify
60 // it under the terms of the GNU General Public License as published by
61 // the Free Software Foundation; either version 2 of the License, or
62 // (at your option) any later version.
63 //
64 // The algorithms that underlie FastJet have required considerable
65 // development and are described in hep-ph/0512210. If you use
66 // FastJet as part of work towards a scientific publication, please
67 // include a citation to the FastJet paper.
68 //
69 // FastJet is distributed in the hope that it will be useful,
70 // but WITHOUT ANY WARRANTY; without even the implied warranty of
71 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
72 // GNU General Public License for more details.
73 //
74 // You should have received a copy of the GNU General Public License
75 // along with FastJet. If not, see <http://www.gnu.org/licenses/>.
76 //----------------------------------------------------------------------
77 //
78 #ifndef __FJCORE_HH__
79 #define __FJCORE_HH__
80 #define __FJCORE__ // remove all the non-core code (a safekeeper)
81 #define __FJCORE_DROP_CGAL // disable CGAL support
82 #ifndef _INCLUDE_FJCORE_CONFIG_AUTO_H
83 #define _INCLUDE_FJCORE_CONFIG_AUTO_H 1
84 #ifndef FJCORE_HAVE_DEMANGLING_SUPPORT
85 #endif
86 #ifndef FJCORE_HAVE_DLFCN_H
87 # define FJCORE_HAVE_DLFCN_H 1
88 #endif
89 #ifndef FJCORE_HAVE_EXECINFO_H
90 #endif
91 #ifndef FJCORE_HAVE_GNUCXX_DEPRECATED
92 #endif
93 #ifndef FJCORE_HAVE_INTTYPES_H
94 # define FJCORE_HAVE_INTTYPES_H 1
95 #endif
96 #ifndef FJCORE_HAVE_LIBM
97 # define FJCORE_HAVE_LIBM 1
98 #endif
99 #ifndef FJCORE_HAVE_MEMORY_H
100 # define FJCORE_HAVE_MEMORY_H 1
101 #endif
102 #ifndef FJCORE_HAVE_STDINT_H
103 # define FJCORE_HAVE_STDINT_H 1
104 #endif
105 #ifndef FJCORE_HAVE_STDLIB_H
106 # define FJCORE_HAVE_STDLIB_H 1
107 #endif
108 #ifndef FJCORE_HAVE_STRINGS_H
109 # define FJCORE_HAVE_STRINGS_H 1
110 #endif
111 #ifndef FJCORE_HAVE_STRING_H
112 # define FJCORE_HAVE_STRING_H 1
113 #endif
114 #ifndef FJCORE_HAVE_SYS_STAT_H
115 # define FJCORE_HAVE_SYS_STAT_H 1
116 #endif
117 #ifndef FJCORE_HAVE_SYS_TYPES_H
118 # define FJCORE_HAVE_SYS_TYPES_H 1
119 #endif
120 #ifndef FJCORE_HAVE_UNISTD_H
121 # define FJCORE_HAVE_UNISTD_H 1
122 #endif
123 #ifndef FJCORE_LT_OBJDIR
124 # define FJCORE_LT_OBJDIR ".libs/"
125 #endif
126 #ifndef FJCORE_PACKAGE
127 # define FJCORE_PACKAGE "fastjet"
128 #endif
129 #ifndef FJCORE_PACKAGE_BUGREPORT
130 # define FJCORE_PACKAGE_BUGREPORT ""
131 #endif
132 #ifndef FJCORE_PACKAGE_NAME
133 # define FJCORE_PACKAGE_NAME "FastJet"
134 #endif
135 #ifndef FJCORE_PACKAGE_STRING
136 # define FJCORE_PACKAGE_STRING "FastJet 3.2.1"
137 #endif
138 #ifndef FJCORE_PACKAGE_TARNAME
139 # define FJCORE_PACKAGE_TARNAME "fastjet"
140 #endif
141 #ifndef FJCORE_PACKAGE_URL
142 # define FJCORE_PACKAGE_URL ""
143 #endif
144 #ifndef FJCORE_PACKAGE_VERSION
145 # define FJCORE_PACKAGE_VERSION "3.2.1"
146 #endif
147 #ifndef FJCORE_STDC_HEADERS
148 # define FJCORE_STDC_HEADERS 1
149 #endif
150 #ifndef FJCORE_VERSION
151 # define FJCORE_VERSION "3.2.1"
152 #endif
153 #ifndef FJCORE_VERSION_MAJOR
154 # define FJCORE_VERSION_MAJOR 3
155 #endif
156 #ifndef FJCORE_VERSION_MINOR
157 # define FJCORE_VERSION_MINOR 2
158 #endif
159 #ifndef FJCORE_VERSION_NUMBER
160 # define FJCORE_VERSION_NUMBER 30201
161 #endif
162 #ifndef FJCORE_VERSION_PATCHLEVEL
163 # define FJCORE_VERSION_PATCHLEVEL 1
164 #endif
165 #endif
166 #ifndef __FJCORE_CONFIG_H__
167 #define __FJCORE_CONFIG_H__
168 #endif // __FJCORE_CONFIG_H__
169 #ifndef __FJCORE_FASTJET_BASE_HH__
170 #define __FJCORE_FASTJET_BASE_HH__
171 // TS : enclose the fjcore namespace inside the Pythia8 one
172 //#define FJCORE_BEGIN_NAMESPACE namespace fjcore {
173 //#define FJCORE_END_NAMESPACE }
174 #define FJCORE_BEGIN_NAMESPACE namespace Pythia8 { namespace fjcore {
175 #define FJCORE_END_NAMESPACE }}
176 #ifdef FJCORE_HAVE_OVERRIDE
177 # define FJCORE_OVERRIDE override
178 #else
179 # define FJCORE_OVERRIDE
180 #endif
181 #endif // __FJCORE_FASTJET_BASE_HH__
182 #ifndef __FJCORE_NUMCONSTS__
183 #define __FJCORE_NUMCONSTS__
184 FJCORE_BEGIN_NAMESPACE // defined in fastjet/internal/base.hh
185 const double pi = 3.141592653589793238462643383279502884197;
186 const double twopi = 6.283185307179586476925286766559005768394;
187 const double pisq = 9.869604401089358618834490999876151135314;
188 const double zeta2 = 1.644934066848226436472415166646025189219;
189 const double zeta3 = 1.202056903159594285399738161511449990765;
190 const double eulergamma = 0.577215664901532860606512090082402431042;
191 const double ln2 = 0.693147180559945309417232121458176568076;
192 FJCORE_END_NAMESPACE
193 #endif // __FJCORE_NUMCONSTS__
194 #ifndef __FJCORE_INTERNAL_IS_BASE_HH__
195 #define __FJCORE_INTERNAL_IS_BASE_HH__
196 FJCORE_BEGIN_NAMESPACE
197 template<typename T, T _t>
198 struct integral_type{
199  static const T value = _t;
200  typedef T value_type;
201  typedef integral_type<T,_t> type;
202 };
203 template<typename T, T _t>
204 const T integral_type<T, _t>::value;
205 typedef integral_type<bool, true> true_type;
206 typedef integral_type<bool, false> false_type;
207 typedef char (&__yes_type)[1]; //< the yes type
208 typedef char (&__no_type) [2]; //< the no type
209 template<typename B, typename D>
210 struct __inheritance_helper{
211 #if !((_MSC_VER !=0 ) && (_MSC_VER == 1310)) // MSVC 7.1
212  template <typename T>
213  static __yes_type check_sig(D const volatile *, T);
214 #else
215  static __yes_type check_sig(D const volatile *, long);
216 #endif
217  static __no_type check_sig(B const volatile *, int);
218 };
219 template<typename B, typename D>
220 struct IsBaseAndDerived{
221 #if ((_MSC_FULL_VER != 0) && (_MSC_FULL_VER >= 140050000))
222 #pragma warning(push)
223 #pragma warning(disable:6334)
224 #endif
225  struct Host{
226 #if !((_MSC_VER !=0 ) && (_MSC_VER == 1310))
227  operator B const volatile *() const;
228 #else
229  operator B const volatile * const&() const;
230 #endif
231  operator D const volatile *();
232  };
233  static const bool value = ((sizeof(B)!=0) &&
234  (sizeof(D)!=0) &&
235  (sizeof(__inheritance_helper<B,D>::check_sig(Host(), 0)) == sizeof(__yes_type)));
236 #if ((_MSC_FULL_VER != 0) && (_MSC_FULL_VER >= 140050000))
237 #pragma warning(pop)
238 #endif
239 };
240 template<class B, class D>
241 B* cast_if_derived(D* d){
242  return IsBaseAndDerived<B,D>::value ? (B*)(d) : 0;
243 }
244 FJCORE_END_NAMESPACE
245 #endif // __IS_BASE_OF_HH__
246 #ifndef __FJCORE_FJCORE_DEPRECATED_HH__
247 #define __FJCORE_FJCORE_DEPRECATED_HH__
248 #if defined(FJCORE_HAVE_CXX14_DEPRECATED)
249 # define FJCORE_DEPRECATED [[deprecated]]
250 # define FJCORE_DEPRECATED_MSG(message) [[deprecated(message)]]
251 #elif defined(FJCORE_HAVE_GNUCXX_DEPRECATED)
252 # define FJCORE_DEPRECATED __attribute__((__deprecated__))
253 # define FJCORE_DEPRECATED_MSG(message) __attribute__((__deprecated__))
254 #else
255 # define FJCORE_DEPRECATED
256 # define FJCORE_DEPRECATED_MSG(message)
257 #endif
258 #endif // __FJCORE_FJCORE_DEPRECATED_HH__
259 #ifndef __FJCORE_SHARED_PTR_HH__
260 #define __FJCORE_SHARED_PTR_HH__
261 #include <cstdlib> // for NULL!!!
262 #ifdef __FJCORE_USETR1SHAREDPTR
263 #include <tr1/memory>
264 #endif // __FJCORE_USETR1SHAREDPTR
265 FJCORE_BEGIN_NAMESPACE // defined in fastjet/internal/base.hh
266 #ifdef __FJCORE_USETR1SHAREDPTR
267 template<class T>
268 class SharedPtr : public std::tr1::shared_ptr<T> {
269 public:
270  SharedPtr() : std::tr1::shared_ptr<T>() {}
271  SharedPtr(T * t) : std::tr1::shared_ptr<T>(t) {}
272  SharedPtr(const SharedPtr<T> & t) : std::tr1::shared_ptr<T>(t) {}
273  #ifdef FJCORE_HAVE_EXPLICIT_FOR_OPERATORS
274  explicit
275  #endif
276  inline operator bool() const {return (this->get()!=NULL);}
277  T* operator ()() const{
278  return this->get(); // automatically returns NULL when out-of-scope
279  }
280 };
281 #else // __FJCORE_USETR1SHAREDPTR
282 template<class T>
283 class SharedPtr{
284 public:
285  class __SharedCountingPtr;
286  SharedPtr() : _ptr(NULL){}
287  template<class Y> explicit SharedPtr(Y* ptr){
288  _ptr = new __SharedCountingPtr(ptr);
289  }
290  SharedPtr(SharedPtr const & share) : _ptr(share._get_container()){
291  if (_ptr!=NULL) ++(*_ptr);
292  }
293  ~SharedPtr(){
294  if (_ptr==NULL) return;
295  _decrease_count();
296  }
297  void reset(){
298  SharedPtr().swap(*this);
299  }
300  template<class Y> void reset(Y * ptr){
301  SharedPtr(ptr).swap(*this);
302  }
303  template<class Y> void reset(SharedPtr<Y> const & share){
304  if (_ptr!=NULL){
305  if (_ptr == share._get_container()) return;
306  _decrease_count();
307  }
308  _ptr = share._get_container(); // Note: automatically set it to NULL if share is empty
309  if (_ptr!=NULL) ++(*_ptr);
310  }
311  SharedPtr& operator=(SharedPtr const & share){
312  reset(share);
313  return *this;
314  }
315  template<class Y> SharedPtr& operator=(SharedPtr<Y> const & share){
316  reset(share);
317  return *this;
318  }
319  FJCORE_DEPRECATED_MSG("Use SharedPtr<T>::get() instead")
320  T* operator ()() const{
321  if (_ptr==NULL) return NULL;
322  return _ptr->get(); // automatically returns NULL when out-of-scope
323  }
324  inline T& operator*() const{
325  return *(_ptr->get());
326  }
327  inline T* operator->() const{
328  if (_ptr==NULL) return NULL;
329  return _ptr->get();
330  }
331  inline T* get() const{
332  if (_ptr==NULL) return NULL;
333  return _ptr->get();
334  }
335  inline bool unique() const{
336  return (use_count()==1);
337  }
338  inline long use_count() const{
339  if (_ptr==NULL) return 0;
340  return _ptr->use_count(); // automatically returns NULL when out-of-scope
341  }
342  #ifdef FJCORE_HAVE_EXPLICIT_FOR_OPERATORS
343  explicit
344  #endif
345  inline operator bool() const{
346  return (get()!=NULL);
347  }
348  inline void swap(SharedPtr & share){
349  __SharedCountingPtr* share_container = share._ptr;
350  share._ptr = _ptr;
351  _ptr = share_container;
352  }
353  void set_count(const long & count){
354  if (_ptr==NULL) return;
355  _ptr->set_count(count);
356  }
357  class __SharedCountingPtr{
358  public:
359  __SharedCountingPtr() : _ptr(NULL), _count(0){}
360  template<class Y> explicit __SharedCountingPtr(Y* ptr) : _ptr(ptr), _count(1){}
361  ~__SharedCountingPtr(){
362  if (_ptr!=NULL){ delete _ptr;}
363  }
364  inline T* get() const {return _ptr;}
365  inline long use_count() const {return _count;}
366  inline long operator++(){return ++_count;}
367  inline long operator--(){return --_count;}
368  inline long operator++(int){return _count++;}
369  inline long operator--(int){return _count--;}
370  void set_count(const long & count){
371  _count = count;
372  }
373  private:
374  T *_ptr;
375  long _count;
376  };
377 private:
378  inline __SharedCountingPtr* _get_container() const{
379  return _ptr;
380  }
381  void _decrease_count(){
382  (*_ptr)--;
383  if (_ptr->use_count()==0)
384  delete _ptr; // that automatically deletes the object itself
385  }
386  __SharedCountingPtr *_ptr;
387 };
388 template<class T,class U>
389 inline bool operator==(SharedPtr<T> const & t, SharedPtr<U> const & u){
390  return t.get() == u.get();
391 }
392 template<class T,class U>
393 inline bool operator!=(SharedPtr<T> const & t, SharedPtr<U> const & u){
394  return t.get() != u.get();
395 }
396 template<class T,class U>
397 inline bool operator<(SharedPtr<T> const & t, SharedPtr<U> const & u){
398  return t.get() < u.get();
399 }
400 template<class T>
401 inline void swap(SharedPtr<T> & a, SharedPtr<T> & b){
402  return a.swap(b);
403 }
404 template<class T>
405 inline T* get_pointer(SharedPtr<T> const & t){
406  return t.get();
407 }
408 #endif // __FJCORE_USETR1SHAREDPTR
409 FJCORE_END_NAMESPACE // defined in fastjet/internal/base.hh
410 #endif // __FJCORE_SHARED_PTR_HH__
411 #ifndef __FJCORE_LIMITEDWARNING_HH__
412 #define __FJCORE_LIMITEDWARNING_HH__
413 #include <iostream>
414 #include <string>
415 #include <list>
416 FJCORE_BEGIN_NAMESPACE
417 class LimitedWarning {
418 public:
419  LimitedWarning() : _max_warn(_max_warn_default), _n_warn_so_far(0), _this_warning_summary(0) {}
420  LimitedWarning(int max_warn_in) : _max_warn(max_warn_in), _n_warn_so_far(0), _this_warning_summary(0) {}
421  void warn(const char * warning) {warn(warning, _default_ostr);}
422  void warn(const std::string & warning) {warn(warning.c_str(), _default_ostr);}
423  void warn(const char * warning, std::ostream * ostr);
424  void warn(const std::string & warning, std::ostream * ostr) {warn(warning.c_str(), ostr);}
425  static void set_default_stream(std::ostream * ostr) {
426  _default_ostr = ostr;
427  }
428  static void set_default_max_warn(int max_warn) {
429  _max_warn_default = max_warn;
430  }
431  int max_warn() const {return _max_warn;}
432  int n_warn_so_far() const {return _n_warn_so_far;}
433  static std::string summary();
434 private:
435  int _max_warn, _n_warn_so_far;
436  static int _max_warn_default;
437  static std::ostream * _default_ostr;
438  typedef std::pair<std::string, unsigned int> Summary;
439  static std::list< Summary > _global_warnings_summary;
440  Summary * _this_warning_summary;
441 };
442 FJCORE_END_NAMESPACE
443 #endif // __FJCORE_LIMITEDWARNING_HH__
444 #ifndef __FJCORE_ERROR_HH__
445 #define __FJCORE_ERROR_HH__
446 #include<iostream>
447 #include<string>
448 #if (!defined(FJCORE_HAVE_EXECINFO_H)) || defined(__FJCORE__)
449 #endif
450 FJCORE_BEGIN_NAMESPACE // defined in fastjet/internal/base.hh
451 class Error {
452 public:
453  Error() {}
454  Error(const std::string & message);
455  virtual ~Error() {}
456  std::string message() const {return _message;}
457  static void set_print_errors(bool print_errors) {_print_errors = print_errors;}
458  static void set_print_backtrace(bool enabled);
459  static void set_default_stream(std::ostream * ostr) {
460  _default_ostr = ostr;
461  }
462 private:
463  std::string _message;
464  static bool _print_errors;
465  static bool _print_backtrace;
466  static std::ostream * _default_ostr;
467 #if (!defined(FJCORE_HAVE_EXECINFO_H)) || defined(__FJCORE__)
468  static LimitedWarning _execinfo_undefined;
469 #endif
470 };
471 class InternalError : public Error{
472 public:
473  InternalError(const std::string & message_in) : Error(std::string("*** CRITICAL INTERNAL FASTJET ERROR *** CONTACT THE AUTHORS *** ") + message_in){ }
474 };
475 FJCORE_END_NAMESPACE
476 #endif // __FJCORE_ERROR_HH__
477 #ifndef __FJCORE_PSEUDOJET_STRUCTURE_BASE_HH__
478 #define __FJCORE_PSEUDOJET_STRUCTURE_BASE_HH__
479 #include <vector>
480 #include <string>
481 FJCORE_BEGIN_NAMESPACE // defined in fastjet/internal/base.hh
482 class PseudoJet;
483 class ClusterSequence;
484 class PseudoJetStructureBase{
485 public:
486  PseudoJetStructureBase(){};
487  virtual ~PseudoJetStructureBase(){};
488  virtual std::string description() const{ return "PseudoJet with an unknown structure"; }
489  virtual bool has_associated_cluster_sequence() const { return false;}
490  virtual const ClusterSequence* associated_cluster_sequence() const;
491  virtual bool has_valid_cluster_sequence() const {return false;}
492  virtual const ClusterSequence * validated_cs() const;
493  virtual bool has_partner(const PseudoJet &reference, PseudoJet &partner) const;
494  virtual bool has_child(const PseudoJet &reference, PseudoJet &child) const;
495  virtual bool has_parents(const PseudoJet &reference, PseudoJet &parent1, PseudoJet &parent2) const;
496  virtual bool object_in_jet(const PseudoJet &reference, const PseudoJet &jet) const;
497  virtual bool has_constituents() const {return false;}
498  virtual std::vector<PseudoJet> constituents(const PseudoJet &reference) const;
499  virtual bool has_exclusive_subjets() const {return false;}
500  virtual std::vector<PseudoJet> exclusive_subjets(const PseudoJet &reference, const double & dcut) const;
501  virtual int n_exclusive_subjets(const PseudoJet &reference, const double & dcut) const;
502  virtual std::vector<PseudoJet> exclusive_subjets_up_to (const PseudoJet &reference, int nsub) const;
503  virtual double exclusive_subdmerge(const PseudoJet &reference, int nsub) const;
504  virtual double exclusive_subdmerge_max(const PseudoJet &reference, int nsub) const;
505  virtual bool has_pieces(const PseudoJet & /* reference */) const {
506  return false;}
507  virtual std::vector<PseudoJet> pieces(const PseudoJet & /* reference */
508  ) const;
509 };
510 FJCORE_END_NAMESPACE
511 #endif // __FJCORE_PSEUDOJET_STRUCTURE_BASE_HH__
512 #ifndef __FJCORE_PSEUDOJET_HH__
513 #define __FJCORE_PSEUDOJET_HH__
514 #include<valarray>
515 #include<vector>
516 #include<cassert>
517 #include<cmath>
518 #include<iostream>
519 FJCORE_BEGIN_NAMESPACE // defined in fastjet/internal/base.hh
520 const double MaxRap = 1e5;
521 const double pseudojet_invalid_phi = -100.0;
522 const double pseudojet_invalid_rap = -1e200;
523 class PseudoJet {
524  public:
525  PseudoJet() : _px(0), _py(0), _pz(0), _E(0) {_finish_init(); _reset_indices();}
526  PseudoJet(const double px, const double py, const double pz, const double E);
527  template <class L> PseudoJet(const L & some_four_vector);
528  PseudoJet(bool /* dummy */) {}
529  virtual ~PseudoJet(){};
530  inline double E() const {return _E;}
531  inline double e() const {return _E;} // like CLHEP
532  inline double px() const {return _px;}
533  inline double py() const {return _py;}
534  inline double pz() const {return _pz;}
535  inline double phi() const {return phi_02pi();}
536  inline double phi_std() const {
537  _ensure_valid_rap_phi();
538  return _phi > pi ? _phi-twopi : _phi;}
539  inline double phi_02pi() const {
540  _ensure_valid_rap_phi();
541  return _phi;
542  }
543  inline double rap() const {
544  _ensure_valid_rap_phi();
545  return _rap;
546  }
547  inline double rapidity() const {return rap();} // like CLHEP
548  double pseudorapidity() const;
549  double eta() const {return pseudorapidity();}
550  inline double pt2() const {return _kt2;}
551  inline double pt() const {return sqrt(_kt2);}
552  inline double perp2() const {return _kt2;} // like CLHEP
553  inline double perp() const {return sqrt(_kt2);} // like CLHEP
554  inline double kt2() const {return _kt2;} // for bkwds compatibility
555  inline double m2() const {return (_E+_pz)*(_E-_pz)-_kt2;}
556  inline double m() const;
557  inline double mperp2() const {return (_E+_pz)*(_E-_pz);}
558  inline double mperp() const {return sqrt(std::abs(mperp2()));}
559  inline double mt2() const {return (_E+_pz)*(_E-_pz);}
560  inline double mt() const {return sqrt(std::abs(mperp2()));}
561  inline double modp2() const {return _kt2+_pz*_pz;}
562  inline double modp() const {return sqrt(_kt2+_pz*_pz);}
563  inline double Et() const {return (_kt2==0) ? 0.0 : _E/sqrt(1.0+_pz*_pz/_kt2);}
564  inline double Et2() const {return (_kt2==0) ? 0.0 : _E*_E/(1.0+_pz*_pz/_kt2);}
565  double operator () (int i) const ;
566  inline double operator [] (int i) const { return (*this)(i); }; // this too
567  double kt_distance(const PseudoJet & other) const;
568  double plain_distance(const PseudoJet & other) const;
569  inline double squared_distance(const PseudoJet & other) const {
570  return plain_distance(other);}
571  inline double delta_R(const PseudoJet & other) const {
572  return sqrt(squared_distance(other));
573  }
574  double delta_phi_to(const PseudoJet & other) const;
575  inline double beam_distance() const {return _kt2;}
576  std::valarray<double> four_mom() const;
577  enum { X=0, Y=1, Z=2, T=3, NUM_COORDINATES=4, SIZE=NUM_COORDINATES };
578  PseudoJet & boost(const PseudoJet & prest);
579  PseudoJet & unboost(const PseudoJet & prest);
580  void operator*=(double);
581  void operator/=(double);
582  void operator+=(const PseudoJet &);
583  void operator-=(const PseudoJet &);
584  inline void reset(double px, double py, double pz, double E);
585  inline void reset(const PseudoJet & psjet) {
586  (*this) = psjet;
587  }
588  template <class L> inline void reset(const L & some_four_vector) {
589  const PseudoJet * pj = fjcore::cast_if_derived<const PseudoJet>(&some_four_vector);
590  if (pj){
591  (*this) = *pj;
592  } else {
593  reset(some_four_vector[0], some_four_vector[1],
594  some_four_vector[2], some_four_vector[3]);
595  }
596  }
597  inline void reset_PtYPhiM(double pt_in, double y_in, double phi_in, double m_in=0.0) {
598  reset_momentum_PtYPhiM(pt_in, y_in, phi_in, m_in);
599  _reset_indices();
600  }
601  inline void reset_momentum(double px, double py, double pz, double E);
602  inline void reset_momentum(const PseudoJet & pj);
603  void reset_momentum_PtYPhiM(double pt, double y, double phi, double m=0.0);
604  template <class L> inline void reset_momentum(const L & some_four_vector) {
605  reset_momentum(some_four_vector[0], some_four_vector[1],
606  some_four_vector[2], some_four_vector[3]);
607  }
608  void set_cached_rap_phi(double rap, double phi);
609  inline int user_index() const {return _user_index;}
610  inline void set_user_index(const int index) {_user_index = index;}
611  class UserInfoBase{
612  public:
613  UserInfoBase(){};
614  virtual ~UserInfoBase(){};
615  };
616  class InexistentUserInfo : public Error {
617  public:
618  InexistentUserInfo();
619  };
620  void set_user_info(UserInfoBase * user_info_in) {
621  _user_info.reset(user_info_in);
622  }
623  template<class L>
624  const L & user_info() const{
625  if (_user_info.get() == 0) throw InexistentUserInfo();
626  return dynamic_cast<const L &>(* _user_info.get());
627  }
628  bool has_user_info() const{
629  return _user_info.get();
630  }
631  template<class L>
632  bool has_user_info() const{
633  return _user_info.get() && dynamic_cast<const L *>(_user_info.get());
634  }
635  const UserInfoBase * user_info_ptr() const{
636  return _user_info.get();
637  }
638  const SharedPtr<UserInfoBase> & user_info_shared_ptr() const{
639  return _user_info;
640  }
641  SharedPtr<UserInfoBase> & user_info_shared_ptr(){
642  return _user_info;
643  }
644  std::string description() const;
645  bool has_associated_cluster_sequence() const;
646  bool has_associated_cs() const {return has_associated_cluster_sequence();}
647  bool has_valid_cluster_sequence() const;
648  bool has_valid_cs() const {return has_valid_cluster_sequence();}
649  const ClusterSequence* associated_cluster_sequence() const;
650  const ClusterSequence* associated_cs() const {return associated_cluster_sequence();}
651  inline const ClusterSequence * validated_cluster_sequence() const {
652  return validated_cs();
653  }
654  const ClusterSequence * validated_cs() const;
655  void set_structure_shared_ptr(const SharedPtr<PseudoJetStructureBase> &structure);
656  bool has_structure() const;
657  const PseudoJetStructureBase* structure_ptr() const;
658  PseudoJetStructureBase* structure_non_const_ptr();
659  const PseudoJetStructureBase* validated_structure_ptr() const;
660  const SharedPtr<PseudoJetStructureBase> & structure_shared_ptr() const;
661  template<typename StructureType>
662  const StructureType & structure() const;
663  template<typename TransformerType>
664  bool has_structure_of() const;
665  template<typename TransformerType>
666  const typename TransformerType::StructureType & structure_of() const;
667  virtual bool has_partner(PseudoJet &partner) const;
668  virtual bool has_child(PseudoJet &child) const;
669  virtual bool has_parents(PseudoJet &parent1, PseudoJet &parent2) const;
670  virtual bool contains(const PseudoJet &constituent) const;
671  virtual bool is_inside(const PseudoJet &jet) const;
672  virtual bool has_constituents() const;
673  virtual std::vector<PseudoJet> constituents() const;
674  virtual bool has_exclusive_subjets() const;
675  std::vector<PseudoJet> exclusive_subjets (const double dcut) const;
676  int n_exclusive_subjets(const double dcut) const;
677  std::vector<PseudoJet> exclusive_subjets (int nsub) const;
678  std::vector<PseudoJet> exclusive_subjets_up_to (int nsub) const;
679  double exclusive_subdmerge(int nsub) const;
680  double exclusive_subdmerge_max(int nsub) const;
681  virtual bool has_pieces() const;
682  virtual std::vector<PseudoJet> pieces() const;
683  inline int cluster_hist_index() const {return _cluster_hist_index;}
684  inline void set_cluster_hist_index(const int index) {_cluster_hist_index = index;}
685  inline int cluster_sequence_history_index() const {
686  return cluster_hist_index();}
687  inline void set_cluster_sequence_history_index(const int index) {
688  set_cluster_hist_index(index);}
689  protected:
690  SharedPtr<PseudoJetStructureBase> _structure;
691  SharedPtr<UserInfoBase> _user_info;
692  private:
693  double _px,_py,_pz,_E;
694  mutable double _phi, _rap;
695  double _kt2;
696  int _cluster_hist_index, _user_index;
697  void _finish_init();
698  void _reset_indices();
699  inline void _ensure_valid_rap_phi() const {
700  if (_phi == pseudojet_invalid_phi) _set_rap_phi();
701  }
702  void _set_rap_phi() const;
703  friend PseudoJet operator*(double, const PseudoJet &);
704 };
705 PseudoJet operator+(const PseudoJet &, const PseudoJet &);
706 PseudoJet operator-(const PseudoJet &, const PseudoJet &);
707 PseudoJet operator*(double, const PseudoJet &);
708 PseudoJet operator*(const PseudoJet &, double);
709 PseudoJet operator/(const PseudoJet &, double);
710 bool operator==(const PseudoJet &, const PseudoJet &);
711 inline bool operator!=(const PseudoJet & a, const PseudoJet & b) {return !(a==b);}
712 bool operator==(const PseudoJet & jet, const double val);
713 inline bool operator==(const double val, const PseudoJet & jet) {return jet == val;}
714 inline bool operator!=(const PseudoJet & a, const double val) {return !(a==val);}
715 inline bool operator!=( const double val, const PseudoJet & a) {return !(a==val);}
716 inline double dot_product(const PseudoJet & a, const PseudoJet & b) {
717  return a.E()*b.E() - a.px()*b.px() - a.py()*b.py() - a.pz()*b.pz();
718 }
719 bool have_same_momentum(const PseudoJet &, const PseudoJet &);
720 PseudoJet PtYPhiM(double pt, double y, double phi, double m = 0.0);
721 std::vector<PseudoJet> sorted_by_pt(const std::vector<PseudoJet> & jets);
722 std::vector<PseudoJet> sorted_by_rapidity(const std::vector<PseudoJet> & jets);
723 std::vector<PseudoJet> sorted_by_E(const std::vector<PseudoJet> & jets);
724 std::vector<PseudoJet> sorted_by_pz(const std::vector<PseudoJet> & jets);
725 void sort_indices(std::vector<int> & indices,
726  const std::vector<double> & values);
727 template<class T> std::vector<T> objects_sorted_by_values(const std::vector<T> & objects,
728  const std::vector<double> & values) {
729  if (objects.size() != values.size()){
730  throw Error("fjcore::objects_sorted_by_values(...): the size of the 'objects' vector must match the size of the 'values' vector");
731  }
732  std::vector<int> indices(values.size());
733  for (size_t i = 0; i < indices.size(); i++) {indices[i] = i;}
734  sort_indices(indices, values);
735  std::vector<T> objects_sorted(objects.size());
736  for (size_t i = 0; i < indices.size(); i++) {
737  objects_sorted[i] = objects[indices[i]];
738  }
739  return objects_sorted;
740 }
741 class IndexedSortHelper {
742 public:
743  inline IndexedSortHelper (const std::vector<double> * reference_values) {
744  _ref_values = reference_values;
745  };
746  inline int operator() (const int i1, const int i2) const {
747  return (*_ref_values)[i1] < (*_ref_values)[i2];
748  };
749 private:
750  const std::vector<double> * _ref_values;
751 };
752 template <class L> inline PseudoJet::PseudoJet(const L & some_four_vector) {
753  reset(some_four_vector);
754 }
755 inline void PseudoJet::_reset_indices() {
756  set_cluster_hist_index(-1);
757  set_user_index(-1);
758  _structure.reset();
759  _user_info.reset();
760 }
761 inline double PseudoJet::m() const {
762  double mm = m2();
763  return mm < 0.0 ? -std::sqrt(-mm) : std::sqrt(mm);
764 }
765 inline void PseudoJet::reset(double px_in, double py_in, double pz_in, double E_in) {
766  _px = px_in;
767  _py = py_in;
768  _pz = pz_in;
769  _E = E_in;
770  _finish_init();
771  _reset_indices();
772 }
773 inline void PseudoJet::reset_momentum(double px_in, double py_in, double pz_in, double E_in) {
774  _px = px_in;
775  _py = py_in;
776  _pz = pz_in;
777  _E = E_in;
778  _finish_init();
779 }
780 inline void PseudoJet::reset_momentum(const PseudoJet & pj) {
781  _px = pj._px ;
782  _py = pj._py ;
783  _pz = pj._pz ;
784  _E = pj._E ;
785  _phi = pj._phi;
786  _rap = pj._rap;
787  _kt2 = pj._kt2;
788 }
789 template<typename StructureType>
790 const StructureType & PseudoJet::structure() const{
791  return dynamic_cast<const StructureType &>(* validated_structure_ptr());
792 }
793 template<typename TransformerType>
794 bool PseudoJet::has_structure_of() const{
795  if (!_structure) return false;
796  return dynamic_cast<const typename TransformerType::StructureType *>(_structure.get()) != 0;
797 }
798 template<typename TransformerType>
799 const typename TransformerType::StructureType & PseudoJet::structure_of() const{
800  if (!_structure)
801  throw Error("Trying to access the structure of a PseudoJet without an associated structure");
802  return dynamic_cast<const typename TransformerType::StructureType &>(*_structure);
803 }
804 PseudoJet join(const std::vector<PseudoJet> & pieces);
805 PseudoJet join(const PseudoJet & j1);
806 PseudoJet join(const PseudoJet & j1, const PseudoJet & j2);
807 PseudoJet join(const PseudoJet & j1, const PseudoJet & j2, const PseudoJet & j3);
808 PseudoJet join(const PseudoJet & j1, const PseudoJet & j2, const PseudoJet & j3, const PseudoJet & j4);
809 FJCORE_END_NAMESPACE
810 #endif // __FJCORE_PSEUDOJET_HH__
811 #ifndef __FJCORE_FUNCTION_OF_PSEUDOJET_HH__
812 #define __FJCORE_FUNCTION_OF_PSEUDOJET_HH__
813 FJCORE_BEGIN_NAMESPACE
814 template<typename TOut>
815 class FunctionOfPseudoJet{
816 public:
817  FunctionOfPseudoJet(){}
818  virtual ~FunctionOfPseudoJet(){}
819  virtual std::string description() const{ return "";}
820  virtual TOut result(const PseudoJet &pj) const = 0;
821  TOut operator()(const PseudoJet &pj) const { return result(pj);}
822  std::vector<TOut> operator()(const std::vector<PseudoJet> &pjs) const {
823  std::vector<TOut> res(pjs.size());
824  for (unsigned int i=0; i<pjs.size(); i++)
825  res[i] = result(pjs[i]);
826  return res;
827  }
828 };
829 FJCORE_END_NAMESPACE
830 #endif // __FJCORE_FUNCTION_OF_PSEUDOJET_HH__
831 #ifndef __FJCORE_SELECTOR_HH__
832 #define __FJCORE_SELECTOR_HH__
833 #include <limits>
834 #include <cmath>
835 FJCORE_BEGIN_NAMESPACE // defined in fastjet/internal/base.hh
836 class Selector;
837 class SelectorWorker {
838 public:
839  virtual ~SelectorWorker() {}
840  virtual bool pass(const PseudoJet & jet) const = 0;
841  virtual void terminator(std::vector<const PseudoJet *> & jets) const {
842  for (unsigned i = 0; i < jets.size(); i++) {
843  if (jets[i] && !pass(*jets[i])) jets[i] = NULL;
844  }
845  }
846  virtual bool applies_jet_by_jet() const {return true;}
847  virtual std::string description() const {return "missing description";}
848  virtual bool takes_reference() const { return false;}
849  virtual void set_reference(const PseudoJet & /*reference*/){
850  throw Error("set_reference(...) cannot be used for a selector worker that does not take a reference");
851  }
852  virtual SelectorWorker* copy(){
853  throw Error("this SelectorWorker has nothing to copy");
854  }
855  virtual void get_rapidity_extent(double & rapmin, double & rapmax) const {
856  rapmax = std::numeric_limits<double>::infinity();
857  rapmin = -rapmax;
858  }
859  virtual bool is_geometric() const { return false;}
860  virtual bool has_finite_area() const;
861  virtual bool has_known_area() const { return false;}
862  virtual double known_area() const{
863  throw Error("this selector has no computable area");
864  }
865 };
866 class Selector{
867 public:
868  Selector() {}
869  Selector(SelectorWorker * worker_in) {_worker.reset(worker_in);}
870  virtual ~Selector(){}
871  bool pass(const PseudoJet & jet) const {
872  if (!validated_worker()->applies_jet_by_jet()) {
873  throw Error("Cannot apply this selector to an individual jet");
874  }
875  return _worker->pass(jet);
876  }
877  bool operator()(const PseudoJet & jet) const {
878  return pass(jet);
879  }
880  unsigned int count(const std::vector<PseudoJet> & jets) const;
881  PseudoJet sum(const std::vector<PseudoJet> & jets) const;
882  double scalar_pt_sum(const std::vector<PseudoJet> & jets) const;
883  void sift(const std::vector<PseudoJet> & jets,
884  std::vector<PseudoJet> & jets_that_pass,
885  std::vector<PseudoJet> & jets_that_fail) const;
886  bool applies_jet_by_jet() const {
887  return validated_worker()->applies_jet_by_jet();
888  }
889  std::vector<PseudoJet> operator()(const std::vector<PseudoJet> & jets) const;
890  virtual void nullify_non_selected(std::vector<const PseudoJet *> & jets) const {
891  validated_worker()->terminator(jets);
892  }
893  void get_rapidity_extent(double &rapmin, double &rapmax) const {
894  return validated_worker()->get_rapidity_extent(rapmin, rapmax);
895  }
896  std::string description() const {
897  return validated_worker()->description();
898  }
899  bool is_geometric() const{
900  return validated_worker()->is_geometric();
901  }
902  bool has_finite_area() const{
903  return validated_worker()->has_finite_area();
904  }
905  const SharedPtr<SelectorWorker> & worker() const {return _worker;}
906  const SelectorWorker* validated_worker() const {
907  const SelectorWorker* worker_ptr = _worker.get();
908  if (worker_ptr == 0) throw InvalidWorker();
909  return worker_ptr;
910  }
911  bool takes_reference() const {
912  return validated_worker()->takes_reference();
913  }
914  const Selector & set_reference(const PseudoJet &reference){
915  if (! validated_worker()->takes_reference()){
916  return *this;
917  }
918  _copy_worker_if_needed();
919  _worker->set_reference(reference);
920  return *this;
921  }
922  class InvalidWorker : public Error {
923  public:
924  InvalidWorker() : Error("Attempt to use Selector with no valid underlying worker") {}
925  };
926  class InvalidArea : public Error {
927  public:
928  InvalidArea() : Error("Attempt to obtain area from Selector for which this is not meaningful") {}
929  };
930  Selector & operator &=(const Selector & b);
931  Selector & operator |=(const Selector & b);
932 protected:
933  void _copy_worker_if_needed(){
934  if (_worker.unique()) return;
935  _worker.reset(_worker->copy());
936  }
937 private:
938  SharedPtr<SelectorWorker> _worker;
939 };
940 Selector SelectorIdentity();
941 Selector operator!(const Selector & s);
942 Selector operator ||(const Selector & s1, const Selector & s2);
943 Selector operator&&(const Selector & s1, const Selector & s2);
944 Selector operator*(const Selector & s1, const Selector & s2);
945 Selector SelectorPtMin(double ptmin);
946 Selector SelectorPtMax(double ptmax);
947 Selector SelectorPtRange(double ptmin, double ptmax);
948 Selector SelectorEtMin(double Etmin);
949 Selector SelectorEtMax(double Etmax);
950 Selector SelectorEtRange(double Etmin, double Etmax);
951 Selector SelectorEMin(double Emin);
952 Selector SelectorEMax(double Emax);
953 Selector SelectorERange(double Emin, double Emax);
954 Selector SelectorMassMin(double Mmin);
955 Selector SelectorMassMax(double Mmax);
956 Selector SelectorMassRange(double Mmin, double Mmax);
957 Selector SelectorRapMin(double rapmin);
958 Selector SelectorRapMax(double rapmax);
959 Selector SelectorRapRange(double rapmin, double rapmax);
960 Selector SelectorAbsRapMin(double absrapmin);
961 Selector SelectorAbsRapMax(double absrapmax);
962 Selector SelectorAbsRapRange(double absrapmin, double absrapmax);
963 Selector SelectorEtaMin(double etamin);
964 Selector SelectorEtaMax(double etamax);
965 Selector SelectorEtaRange(double etamin, double etamax);
966 Selector SelectorAbsEtaMin(double absetamin);
967 Selector SelectorAbsEtaMax(double absetamax);
968 Selector SelectorAbsEtaRange(double absetamin, double absetamax);
969 Selector SelectorPhiRange(double phimin, double phimax);
970 Selector SelectorRapPhiRange(double rapmin, double rapmax, double phimin, double phimax);
971 Selector SelectorNHardest(unsigned int n);
972 Selector SelectorCircle(const double radius);
973 Selector SelectorDoughnut(const double radius_in, const double radius_out);
974 Selector SelectorStrip(const double half_width);
975 Selector SelectorRectangle(const double half_rap_width, const double half_phi_width);
976 Selector SelectorPtFractionMin(double fraction);
977 Selector SelectorIsZero();
978 FJCORE_END_NAMESPACE // defined in fastjet/internal/base.hh
979 #endif // __FJCORE_SELECTOR_HH__
980 #ifndef __FJCORE_JETDEFINITION_HH__
981 #define __FJCORE_JETDEFINITION_HH__
982 #include<cassert>
983 #include<string>
984 #include<memory>
985 FJCORE_BEGIN_NAMESPACE // defined in fastjet/internal/base.hh
986 std::string fastjet_version_string();
987 enum Strategy {
988  N2MHTLazy9AntiKtSeparateGhosts = -10,
989  N2MHTLazy9 = -7,
990  N2MHTLazy25 = -6,
991  N2MHTLazy9Alt = -5,
992  N2MinHeapTiled = -4,
993  N2Tiled = -3,
994  N2PoorTiled = -2,
995  N2Plain = -1,
996  N3Dumb = 0,
997  Best = 1,
998  NlnN = 2,
999  NlnN3pi = 3,
1000  NlnN4pi = 4,
1001  NlnNCam4pi = 14,
1002  NlnNCam2pi2R = 13,
1003  NlnNCam = 12, // 2piMultD
1004  BestFJ30 = 21,
1005  plugin_strategy = 999
1006 };
1007 enum JetAlgorithm {
1008  kt_algorithm=0,
1009  cambridge_algorithm=1,
1010  antikt_algorithm=2,
1011  genkt_algorithm=3,
1012  cambridge_for_passive_algorithm=11,
1013  genkt_for_passive_algorithm=13,
1014  ee_kt_algorithm=50,
1015  ee_genkt_algorithm=53,
1016  plugin_algorithm = 99,
1017  undefined_jet_algorithm = 999
1018 };
1019 typedef JetAlgorithm JetFinder;
1020 const JetAlgorithm aachen_algorithm = cambridge_algorithm;
1021 const JetAlgorithm cambridge_aachen_algorithm = cambridge_algorithm;
1022 enum RecombinationScheme {
1023  E_scheme=0,
1024  pt_scheme=1,
1025  pt2_scheme=2,
1026  Et_scheme=3,
1027  Et2_scheme=4,
1028  BIpt_scheme=5,
1029  BIpt2_scheme=6,
1030  WTA_pt_scheme=7,
1031  WTA_modp_scheme=8,
1032  external_scheme = 99
1033 };
1034 class ClusterSequence;
1035 class JetDefinition {
1036 public:
1037  class Plugin;
1038  class Recombiner;
1039  JetDefinition(JetAlgorithm jet_algorithm_in,
1040  double R_in,
1041  RecombinationScheme recomb_scheme_in = E_scheme,
1042  Strategy strategy_in = Best) {
1043  *this = JetDefinition(jet_algorithm_in, R_in, recomb_scheme_in, strategy_in, 1);
1044  }
1045  JetDefinition(JetAlgorithm jet_algorithm_in,
1046  RecombinationScheme recomb_scheme_in = E_scheme,
1047  Strategy strategy_in = Best) {
1048  double dummyR = 0.0;
1049  *this = JetDefinition(jet_algorithm_in, dummyR, recomb_scheme_in, strategy_in, 0);
1050  }
1051  JetDefinition(JetAlgorithm jet_algorithm_in,
1052  double R_in,
1053  double xtra_param_in,
1054  RecombinationScheme recomb_scheme_in = E_scheme,
1055  Strategy strategy_in = Best) {
1056  *this = JetDefinition(jet_algorithm_in, R_in, recomb_scheme_in, strategy_in, 2);
1057  set_extra_param(xtra_param_in);
1058  }
1059  JetDefinition(JetAlgorithm jet_algorithm_in,
1060  double R_in,
1061  const Recombiner * recombiner_in,
1062  Strategy strategy_in = Best) {
1063  *this = JetDefinition(jet_algorithm_in, R_in, external_scheme, strategy_in);
1064  _recombiner = recombiner_in;
1065  }
1066  JetDefinition(JetAlgorithm jet_algorithm_in,
1067  const Recombiner * recombiner_in,
1068  Strategy strategy_in = Best) {
1069  *this = JetDefinition(jet_algorithm_in, external_scheme, strategy_in);
1070  _recombiner = recombiner_in;
1071  }
1072  JetDefinition(JetAlgorithm jet_algorithm_in,
1073  double R_in,
1074  double xtra_param_in,
1075  const Recombiner * recombiner_in,
1076  Strategy strategy_in = Best) {
1077  *this = JetDefinition(jet_algorithm_in, R_in, xtra_param_in, external_scheme, strategy_in);
1078  _recombiner = recombiner_in;
1079  }
1080  JetDefinition() {
1081  *this = JetDefinition(undefined_jet_algorithm, 1.0);
1082  }
1083  JetDefinition(const Plugin * plugin_in) {
1084  _plugin = plugin_in;
1085  _strategy = plugin_strategy;
1086  _Rparam = _plugin->R();
1087  _extra_param = 0.0; // a dummy value to keep static code checkers happy
1088  _jet_algorithm = plugin_algorithm;
1089  set_recombination_scheme(E_scheme);
1090  }
1091  JetDefinition(JetAlgorithm jet_algorithm_in,
1092  double R_in,
1093  RecombinationScheme recomb_scheme_in,
1094  Strategy strategy_in,
1095  int nparameters_in);
1096  FJCORE_DEPRECATED_MSG("This argument ordering is deprecated. Use JetDefinition(alg, R, strategy, scheme[, n_parameters]) instead")
1097  JetDefinition(JetAlgorithm jet_algorithm_in,
1098  double R_in,
1099  Strategy strategy_in,
1100  RecombinationScheme recomb_scheme_in = E_scheme,
1101  int nparameters_in = 1){
1102  (*this) = JetDefinition(jet_algorithm_in,R_in,recomb_scheme_in,strategy_in,nparameters_in);
1103  }
1104  template <class L>
1105  std::vector<PseudoJet> operator()(const std::vector<L> & particles) const;
1106  static const double max_allowable_R; //= 1000.0;
1107  void set_recombination_scheme(RecombinationScheme);
1108  void set_recombiner(const Recombiner * recomb) {
1109  if (_shared_recombiner) _shared_recombiner.reset(recomb);
1110  _recombiner = recomb;
1111  _default_recombiner = DefaultRecombiner(external_scheme);
1112  }
1113  void set_recombiner(const JetDefinition &other_jet_def);
1114  void delete_recombiner_when_unused();
1115  const Plugin * plugin() const {return _plugin;};
1116  void delete_plugin_when_unused();
1117  JetAlgorithm jet_algorithm () const {return _jet_algorithm ;}
1118  JetAlgorithm jet_finder () const {return _jet_algorithm ;}
1119  double R () const {return _Rparam ;}
1120  double extra_param () const {return _extra_param ;}
1121  Strategy strategy () const {return _strategy ;}
1122  RecombinationScheme recombination_scheme() const {
1123  return _default_recombiner.scheme();}
1124  void set_jet_algorithm(JetAlgorithm njf) {_jet_algorithm = njf;}
1125  void set_jet_finder(JetAlgorithm njf) {_jet_algorithm = njf;}
1126  void set_extra_param(double xtra_param) {_extra_param = xtra_param;}
1127  const Recombiner * recombiner() const {
1128  return _recombiner == 0 ? & _default_recombiner : _recombiner;}
1129  bool has_same_recombiner(const JetDefinition &other_jd) const;
1130  bool is_spherical() const;
1131  std::string description() const;
1132  std::string description_no_recombiner() const;
1133  static std::string algorithm_description(const JetAlgorithm jet_alg);
1134  static unsigned int n_parameters_for_algorithm(const JetAlgorithm jet_alg);
1135 public:
1136  class Recombiner {
1137  public:
1138  virtual std::string description() const = 0;
1139  virtual void recombine(const PseudoJet & pa, const PseudoJet & pb,
1140  PseudoJet & pab) const = 0;
1141  virtual void preprocess(PseudoJet & ) const {};
1142  virtual ~Recombiner() {};
1143  inline void plus_equal(PseudoJet & pa, const PseudoJet & pb) const {
1144  PseudoJet pres;
1145  recombine(pa,pb,pres);
1146  pa = pres;
1147  }
1148  };
1149  class DefaultRecombiner : public Recombiner {
1150  public:
1151  DefaultRecombiner(RecombinationScheme recomb_scheme = E_scheme) :
1152  _recomb_scheme(recomb_scheme) {}
1153  virtual std::string description() const FJCORE_OVERRIDE;
1154  virtual void recombine(const PseudoJet & pa, const PseudoJet & pb,
1155  PseudoJet & pab) const FJCORE_OVERRIDE;
1156  virtual void preprocess(PseudoJet & p) const FJCORE_OVERRIDE;
1157  RecombinationScheme scheme() const {return _recomb_scheme;}
1158  private:
1159  RecombinationScheme _recomb_scheme;
1160  };
1161  class Plugin{
1162  public:
1163  virtual std::string description() const = 0;
1164  virtual void run_clustering(ClusterSequence &) const = 0;
1165  virtual double R() const = 0;
1166  virtual bool supports_ghosted_passive_areas() const {return false;}
1167  virtual void set_ghost_separation_scale(double scale) const;
1168  virtual double ghost_separation_scale() const {return 0.0;}
1169  virtual bool exclusive_sequence_meaningful() const {return false;}
1170  virtual bool is_spherical() const {return false;}
1171  virtual ~Plugin() {};
1172  };
1173 private:
1174  JetAlgorithm _jet_algorithm;
1175  double _Rparam;
1176  double _extra_param ;
1177  Strategy _strategy ;
1178  const Plugin * _plugin;
1179  SharedPtr<const Plugin> _plugin_shared;
1180  DefaultRecombiner _default_recombiner;
1181  const Recombiner * _recombiner;
1182  SharedPtr<const Recombiner> _shared_recombiner;
1183 };
1184 PseudoJet join(const std::vector<PseudoJet> & pieces, const JetDefinition::Recombiner & recombiner);
1185 PseudoJet join(const PseudoJet & j1,
1186  const JetDefinition::Recombiner & recombiner);
1187 PseudoJet join(const PseudoJet & j1, const PseudoJet & j2,
1188  const JetDefinition::Recombiner & recombiner);
1189 PseudoJet join(const PseudoJet & j1, const PseudoJet & j2, const PseudoJet & j3,
1190  const JetDefinition::Recombiner & recombiner);
1191 PseudoJet join(const PseudoJet & j1, const PseudoJet & j2, const PseudoJet & j3, const PseudoJet & j4,
1192  const JetDefinition::Recombiner & recombiner);
1193 FJCORE_END_NAMESPACE
1194 #endif // __FJCORE_JETDEFINITION_HH__
1195 #ifndef __FJCORE_COMPOSITEJET_STRUCTURE_HH__
1196 #define __FJCORE_COMPOSITEJET_STRUCTURE_HH__
1197 FJCORE_BEGIN_NAMESPACE // defined in fastjet/internal/base.hh
1198 class CompositeJetStructure : public PseudoJetStructureBase{
1199 public:
1200  CompositeJetStructure() : _area_4vector_ptr(0){};
1201  CompositeJetStructure(const std::vector<PseudoJet> & initial_pieces,
1202  const JetDefinition::Recombiner * recombiner = 0);
1203  virtual ~CompositeJetStructure(){
1204  if (_area_4vector_ptr) delete _area_4vector_ptr;
1205  };
1206  virtual std::string description() const FJCORE_OVERRIDE;
1207  virtual bool has_constituents() const FJCORE_OVERRIDE;
1208  virtual std::vector<PseudoJet> constituents(const PseudoJet &jet) const FJCORE_OVERRIDE;
1209  virtual bool has_pieces(const PseudoJet & /*jet*/) const FJCORE_OVERRIDE {return true;}
1210  virtual std::vector<PseudoJet> pieces(const PseudoJet &jet) const FJCORE_OVERRIDE;
1211 protected:
1212  std::vector<PseudoJet> _pieces;
1213  PseudoJet * _area_4vector_ptr;
1214 };
1215 template<typename T> PseudoJet join(const std::vector<PseudoJet> & pieces){
1216  PseudoJet result(0.0,0.0,0.0,0.0);
1217  for (unsigned int i=0; i<pieces.size(); i++){
1218  const PseudoJet it = pieces[i];
1219  result += it;
1220  }
1221  T *cj_struct = new T(pieces);
1222  result.set_structure_shared_ptr(SharedPtr<PseudoJetStructureBase>(cj_struct));
1223  return result;
1224 }
1225 template<typename T> PseudoJet join(const PseudoJet & j1){
1226  return join<T>(std::vector<PseudoJet>(1,j1));
1227 }
1228 template<typename T> PseudoJet join(const PseudoJet & j1, const PseudoJet & j2){
1229  std::vector<PseudoJet> pieces;
1230  pieces.push_back(j1);
1231  pieces.push_back(j2);
1232  return join<T>(pieces);
1233 }
1234 template<typename T> PseudoJet join(const PseudoJet & j1, const PseudoJet & j2,
1235  const PseudoJet & j3){
1236  std::vector<PseudoJet> pieces;
1237  pieces.push_back(j1);
1238  pieces.push_back(j2);
1239  pieces.push_back(j3);
1240  return join<T>(pieces);
1241 }
1242 template<typename T> PseudoJet join(const PseudoJet & j1, const PseudoJet & j2,
1243  const PseudoJet & j3, const PseudoJet & j4){
1244  std::vector<PseudoJet> pieces;
1245  pieces.push_back(j1);
1246  pieces.push_back(j2);
1247  pieces.push_back(j3);
1248  pieces.push_back(j4);
1249  return join<T>(pieces);
1250 }
1251 template<typename T> PseudoJet join(const std::vector<PseudoJet> & pieces,
1252  const JetDefinition::Recombiner & recombiner){
1253  PseudoJet result;
1254  if (pieces.size()>0){
1255  result = pieces[0];
1256  for (unsigned int i=1; i<pieces.size(); i++){
1257  recombiner.plus_equal(result, pieces[i]);
1258  }
1259  }
1260  T *cj_struct = new T(pieces, &recombiner);
1261  result.set_structure_shared_ptr(SharedPtr<PseudoJetStructureBase>(cj_struct));
1262  return result;
1263 }
1264 template<typename T> PseudoJet join(const PseudoJet & j1,
1265  const JetDefinition::Recombiner & recombiner){
1266  return join<T>(std::vector<PseudoJet>(1,j1), recombiner);
1267 }
1268 template<typename T> PseudoJet join(const PseudoJet & j1, const PseudoJet & j2,
1269  const JetDefinition::Recombiner & recombiner){
1270  std::vector<PseudoJet> pieces;
1271  pieces.reserve(2);
1272  pieces.push_back(j1);
1273  pieces.push_back(j2);
1274  return join<T>(pieces, recombiner);
1275 }
1276 template<typename T> PseudoJet join(const PseudoJet & j1, const PseudoJet & j2,
1277  const PseudoJet & j3,
1278  const JetDefinition::Recombiner & recombiner){
1279  std::vector<PseudoJet> pieces;
1280  pieces.reserve(3);
1281  pieces.push_back(j1);
1282  pieces.push_back(j2);
1283  pieces.push_back(j3);
1284  return join<T>(pieces, recombiner);
1285 }
1286 template<typename T> PseudoJet join(const PseudoJet & j1, const PseudoJet & j2,
1287  const PseudoJet & j3, const PseudoJet & j4,
1288  const JetDefinition::Recombiner & recombiner){
1289  std::vector<PseudoJet> pieces;
1290  pieces.reserve(4);
1291  pieces.push_back(j1);
1292  pieces.push_back(j2);
1293  pieces.push_back(j3);
1294  pieces.push_back(j4);
1295  return join<T>(pieces, recombiner);
1296 }
1297 FJCORE_END_NAMESPACE // defined in fastjet/internal/base.hh
1298 #endif // __FJCORE_MERGEDJET_STRUCTURE_HH__
1299 #ifndef __FJCORE_CLUSTER_SEQUENCE_STRUCTURE_HH__
1300 #define __FJCORE_CLUSTER_SEQUENCE_STRUCTURE_HH__
1301 #include <vector>
1302 FJCORE_BEGIN_NAMESPACE // defined in fastjet/internal/base.hh
1303 class ClusterSequenceStructure : public PseudoJetStructureBase{
1304 public:
1305  ClusterSequenceStructure() : _associated_cs(NULL){}
1306  ClusterSequenceStructure(const ClusterSequence *cs){
1307  set_associated_cs(cs);
1308  };
1309  virtual ~ClusterSequenceStructure();
1310  virtual std::string description() const FJCORE_OVERRIDE{
1311  return "PseudoJet with an associated ClusterSequence";
1312  }
1313  virtual bool has_associated_cluster_sequence() const FJCORE_OVERRIDE{ return true;}
1314  virtual const ClusterSequence* associated_cluster_sequence() const FJCORE_OVERRIDE;
1315  virtual bool has_valid_cluster_sequence() const FJCORE_OVERRIDE;
1316  virtual const ClusterSequence * validated_cs() const FJCORE_OVERRIDE;
1317  virtual void set_associated_cs(const ClusterSequence * new_cs){
1318  _associated_cs = new_cs;
1319  }
1320  virtual bool has_partner(const PseudoJet &reference, PseudoJet &partner) const FJCORE_OVERRIDE;
1321  virtual bool has_child(const PseudoJet &reference, PseudoJet &child) const FJCORE_OVERRIDE;
1322  virtual bool has_parents(const PseudoJet &reference, PseudoJet &parent1, PseudoJet &parent2) const FJCORE_OVERRIDE;
1323  virtual bool object_in_jet(const PseudoJet &reference, const PseudoJet &jet) const FJCORE_OVERRIDE;
1324  virtual bool has_constituents() const FJCORE_OVERRIDE;
1325  virtual std::vector<PseudoJet> constituents(const PseudoJet &reference) const FJCORE_OVERRIDE;
1326  virtual bool has_exclusive_subjets() const FJCORE_OVERRIDE;
1327  virtual std::vector<PseudoJet> exclusive_subjets(const PseudoJet &reference, const double & dcut) const FJCORE_OVERRIDE;
1328  virtual int n_exclusive_subjets(const PseudoJet &reference, const double & dcut) const FJCORE_OVERRIDE;
1329  virtual std::vector<PseudoJet> exclusive_subjets_up_to (const PseudoJet &reference, int nsub) const FJCORE_OVERRIDE;
1330  virtual double exclusive_subdmerge(const PseudoJet &reference, int nsub) const FJCORE_OVERRIDE;
1331  virtual double exclusive_subdmerge_max(const PseudoJet &reference, int nsub) const FJCORE_OVERRIDE;
1332  virtual bool has_pieces(const PseudoJet &reference) const FJCORE_OVERRIDE;
1333  virtual std::vector<PseudoJet> pieces(const PseudoJet &reference) const FJCORE_OVERRIDE;
1334 protected:
1335  const ClusterSequence *_associated_cs;
1336 };
1337 FJCORE_END_NAMESPACE
1338 #endif // __FJCORE_CLUSTER_SEQUENCE_STRUCTURE_HH__
1339 #ifndef __FJCORE_CLUSTERSEQUENCE_HH__
1340 #define __FJCORE_CLUSTERSEQUENCE_HH__
1341 #include<vector>
1342 #include<map>
1343 #include<memory>
1344 #include<cassert>
1345 #include<iostream>
1346 #include<string>
1347 #include<set>
1348 #include<cmath> // needed to get double std::abs(double)
1349 FJCORE_BEGIN_NAMESPACE // defined in fastjet/internal/base.hh
1350 class ClusterSequenceStructure;
1351 class DynamicNearestNeighbours;
1352 class ClusterSequence {
1353  public:
1354  ClusterSequence () : _deletes_self_when_unused(false) {}
1355  template<class L> ClusterSequence (
1356  const std::vector<L> & pseudojets,
1357  const JetDefinition & jet_def,
1358  const bool & writeout_combinations = false);
1359  ClusterSequence (const ClusterSequence & cs) : _deletes_self_when_unused(false) {
1360  transfer_from_sequence(cs);
1361  }
1362  ClusterSequence & operator=(const ClusterSequence & cs);
1363  virtual ~ClusterSequence (); //{}
1364  std::vector<PseudoJet> inclusive_jets (const double ptmin = 0.0) const;
1365  int n_exclusive_jets (const double dcut) const;
1366  std::vector<PseudoJet> exclusive_jets (const double dcut) const;
1367  std::vector<PseudoJet> exclusive_jets (const int njets) const;
1368  std::vector<PseudoJet> exclusive_jets_up_to (const int njets) const;
1369  double exclusive_dmerge (const int njets) const;
1370  double exclusive_dmerge_max (const int njets) const;
1371  double exclusive_ymerge (int njets) const {return exclusive_dmerge(njets) / Q2();}
1372  double exclusive_ymerge_max (int njets) const {return exclusive_dmerge_max(njets)/Q2();}
1373  int n_exclusive_jets_ycut (double ycut) const {return n_exclusive_jets(ycut*Q2());}
1374  std::vector<PseudoJet> exclusive_jets_ycut (double ycut) const {
1375  int njets = n_exclusive_jets_ycut(ycut);
1376  return exclusive_jets(njets);
1377  }
1378  std::vector<PseudoJet> exclusive_subjets (const PseudoJet & jet,
1379  const double dcut) const;
1380  int n_exclusive_subjets(const PseudoJet & jet,
1381  const double dcut) const;
1382  std::vector<PseudoJet> exclusive_subjets (const PseudoJet & jet,
1383  int nsub) const;
1384  std::vector<PseudoJet> exclusive_subjets_up_to (const PseudoJet & jet,
1385  int nsub) const;
1386  double exclusive_subdmerge(const PseudoJet & jet, int nsub) const;
1387  double exclusive_subdmerge_max(const PseudoJet & jet, int nsub) const;
1388  double Q() const {return _Qtot;}
1389  double Q2() const {return _Qtot*_Qtot;}
1390  bool object_in_jet(const PseudoJet & object, const PseudoJet & jet) const;
1391  bool has_parents(const PseudoJet & jet, PseudoJet & parent1,
1392  PseudoJet & parent2) const;
1393  bool has_child(const PseudoJet & jet, PseudoJet & child) const;
1394  bool has_child(const PseudoJet & jet, const PseudoJet * & childp) const;
1395  bool has_partner(const PseudoJet & jet, PseudoJet & partner) const;
1396  std::vector<PseudoJet> constituents (const PseudoJet & jet) const;
1397  void print_jets_for_root(const std::vector<PseudoJet> & jets,
1398  std::ostream & ostr = std::cout) const;
1399  void print_jets_for_root(const std::vector<PseudoJet> & jets,
1400  const std::string & filename,
1401  const std::string & comment = "") const;
1402  void add_constituents (const PseudoJet & jet,
1403  std::vector<PseudoJet> & subjet_vector) const;
1404  inline Strategy strategy_used () const {return _strategy;}
1405  std::string strategy_string () const {return strategy_string(_strategy);}
1406  std::string strategy_string (Strategy strategy_in) const;
1407  const JetDefinition & jet_def() const {return _jet_def;}
1408  void delete_self_when_unused();
1409  bool will_delete_self_when_unused() const {return _deletes_self_when_unused;}
1410  void signal_imminent_self_deletion() const;
1411  double jet_scale_for_algorithm(const PseudoJet & jet) const;
1412  void plugin_record_ij_recombination(int jet_i, int jet_j, double dij,
1413  int & newjet_k) {
1414  assert(plugin_activated());
1415  _do_ij_recombination_step(jet_i, jet_j, dij, newjet_k);
1416  }
1417  void plugin_record_ij_recombination(int jet_i, int jet_j, double dij,
1418  const PseudoJet & newjet,
1419  int & newjet_k);
1420  void plugin_record_iB_recombination(int jet_i, double diB) {
1421  assert(plugin_activated());
1422  _do_iB_recombination_step(jet_i, diB);
1423  }
1424  class Extras {
1425  public:
1426  virtual ~Extras() {}
1427  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";}
1428  };
1429  inline void plugin_associate_extras(Extras * extras_in) {
1430  _extras.reset(extras_in);
1431  }
1432 // TS : fully disable use of auto_ptr (already deprecated).
1433 //#ifdef FJCORE_HAVE_AUTO_PTR_INTERFACE
1434 // FJCORE_DEPRECATED_MSG("Please use ClusterSequence::plugin_associate_extras(Extras * extras_in)) instead")
1435 // inline void plugin_associate_extras(std::auto_ptr<Extras> extras_in){
1436 // _extras.reset(extras_in.release());
1437 // }
1438 //#endif
1439  inline bool plugin_activated() const {return _plugin_activated;}
1440  const Extras * extras() const {return _extras.get();}
1441  template<class GBJ> void plugin_simple_N2_cluster () {
1442  assert(plugin_activated());
1443  _simple_N2_cluster<GBJ>();
1444  }
1445 public:
1446  struct history_element{
1447  int parent1;
1448  int parent2;
1449  int child;
1450  int jetp_index;
1454  double dij;
1455  double max_dij_so_far;
1457  };
1459  enum JetType {Invalid=-3, InexistentParent = -2, BeamJet = -1};
1460  const std::vector<PseudoJet> & jets() const;
1461  const std::vector<history_element> & history() const;
1462  unsigned int n_particles() const;
1463  std::vector<int> particle_jet_indices(const std::vector<PseudoJet> &) const;
1464  std::vector<int> unique_history_order() const;
1465  std::vector<PseudoJet> unclustered_particles() const;
1466  std::vector<PseudoJet> childless_pseudojets() const;
1467  bool contains(const PseudoJet & object) const;
1468  void transfer_from_sequence(const ClusterSequence & from_seq,
1469  const FunctionOfPseudoJet<PseudoJet> * action_on_jets = 0);
1470  const SharedPtr<PseudoJetStructureBase> & structure_shared_ptr() const{
1471  return _structure_shared_ptr;
1472  }
1473  typedef ClusterSequenceStructure StructureType;
1474  static void print_banner();
1475  static void set_fastjet_banner_stream(std::ostream * ostr) {_fastjet_banner_ostr = ostr;}
1476  static std::ostream * fastjet_banner_stream() {return _fastjet_banner_ostr;}
1477 private:
1478  static std::ostream * _fastjet_banner_ostr;
1479 protected:
1480  JetDefinition _jet_def;
1481  template<class L> void _transfer_input_jets(
1482  const std::vector<L> & pseudojets);
1483  void _initialise_and_run (const JetDefinition & jet_def,
1484  const bool & writeout_combinations);
1485  void _initialise_and_run_no_decant();
1486  void _decant_options(const JetDefinition & jet_def,
1487  const bool & writeout_combinations);
1488  void _decant_options_partial();
1489  void _fill_initial_history();
1490  void _do_ij_recombination_step(const int jet_i, const int jet_j,
1491  const double dij, int & newjet_k);
1492  void _do_iB_recombination_step(const int jet_i, const double diB);
1493  void _set_structure_shared_ptr(PseudoJet & j);
1494  void _update_structure_use_count();
1495  Strategy _best_strategy() const;
1496  class _Parabola {
1497  public:
1498  _Parabola(double a, double b, double c) : _a(a), _b(b), _c(c) {}
1499  inline double operator()(const double R) const {return _c*(_a*R*R + _b*R + 1);}
1500  private:
1501  double _a, _b, _c;
1502  };
1503  class _Line {
1504  public:
1505  _Line(double a, double b) : _a(a), _b(b) {}
1506  inline double operator()(const double R) const {return _a*R + _b;}
1507  private:
1508  double _a, _b;
1509  };
1510  std::vector<PseudoJet> _jets;
1511  std::vector<history_element> _history;
1512  void get_subhist_set(std::set<const history_element*> & subhist,
1513  const PseudoJet & jet, double dcut, int maxjet) const;
1514  bool _writeout_combinations;
1515  int _initial_n;
1516  double _Rparam, _R2, _invR2;
1517  double _Qtot;
1518  Strategy _strategy;
1519  JetAlgorithm _jet_algorithm;
1520  SharedPtr<PseudoJetStructureBase> _structure_shared_ptr; //< will actually be of type ClusterSequenceStructure
1521  int _structure_use_count_after_construction; //< info of use when CS handles its own memory
1522  mutable bool _deletes_self_when_unused;
1523  private:
1524  bool _plugin_activated;
1525  SharedPtr<Extras> _extras; // things the plugin might want to add
1526  void _really_dumb_cluster ();
1527  void _delaunay_cluster ();
1528  template<class BJ> void _simple_N2_cluster ();
1529  void _tiled_N2_cluster ();
1530  void _faster_tiled_N2_cluster ();
1531  void _minheap_faster_tiled_N2_cluster();
1532  void _CP2DChan_cluster();
1533  void _CP2DChan_cluster_2pi2R ();
1534  void _CP2DChan_cluster_2piMultD ();
1535  void _CP2DChan_limited_cluster(double D);
1536  void _do_Cambridge_inclusive_jets();
1537  void _fast_NsqrtN_cluster();
1538  void _add_step_to_history( //const int step_number,
1539  const int parent1,
1540  const int parent2, const int jetp_index,
1541  const double dij);
1542  void _extract_tree_children(int pos, std::valarray<bool> &,
1543  const std::valarray<int> &, std::vector<int> &) const;
1544  void _extract_tree_parents (int pos, std::valarray<bool> &,
1545  const std::valarray<int> &, std::vector<int> &) const;
1546  typedef std::pair<int,int> TwoVertices;
1547  typedef std::pair<double,TwoVertices> DijEntry;
1548  typedef std::multimap<double,TwoVertices> DistMap;
1549  void _add_ktdistance_to_map(const int ii,
1550  DistMap & DijMap,
1551  const DynamicNearestNeighbours * DNN);
1552  static bool _first_time;
1553  static LimitedWarning _exclusive_warnings;
1554  static LimitedWarning _changed_strategy_warning;
1555  struct BriefJet {
1556  double eta, phi, kt2, NN_dist;
1557  BriefJet * NN;
1558  int _jets_index;
1559  };
1560  class TiledJet {
1561  public:
1562  double eta, phi, kt2, NN_dist;
1563  TiledJet * NN, *previous, * next;
1564  int _jets_index, tile_index, diJ_posn;
1565  inline void label_minheap_update_needed() {diJ_posn = 1;}
1566  inline void label_minheap_update_done() {diJ_posn = 0;}
1567  inline bool minheap_update_needed() const {return diJ_posn==1;}
1568  };
1569  template <class J> void _bj_set_jetinfo( J * const jet,
1570  const int _jets_index) const;
1571  void _bj_remove_from_tiles( TiledJet * const jet) const;
1572  template <class J> double _bj_dist(const J * const jeta,
1573  const J * const jetb) const;
1574  template <class J> double _bj_diJ(const J * const jeta) const;
1575  template <class J> inline J * _bj_of_hindex(
1576  const int hist_index,
1577  J * const head, J * const tail)
1578  const {
1579  J * res;
1580  for(res = head; res<tail; res++) {
1581  if (_jets[res->_jets_index].cluster_hist_index() == hist_index) {break;}
1582  }
1583  return res;
1584  }
1585  template <class J> void _bj_set_NN_nocross(J * const jeta,
1586  J * const head, const J * const tail) const;
1587  template <class J> void _bj_set_NN_crosscheck(J * const jeta,
1588  J * const head, const J * const tail) const;
1589  static const int n_tile_neighbours = 9;
1590  struct Tile {
1591  Tile * begin_tiles[n_tile_neighbours];
1592  Tile ** surrounding_tiles;
1593  Tile ** RH_tiles;
1594  Tile ** end_tiles;
1595  TiledJet * head;
1596  bool tagged;
1597  };
1598  std::vector<Tile> _tiles;
1599  double _tiles_eta_min, _tiles_eta_max;
1600  double _tile_size_eta, _tile_size_phi;
1601  int _n_tiles_phi,_tiles_ieta_min,_tiles_ieta_max;
1602  inline int _tile_index (int ieta, int iphi) const {
1603  return (ieta-_tiles_ieta_min)*_n_tiles_phi
1604  + (iphi+_n_tiles_phi) % _n_tiles_phi;
1605  }
1606  int _tile_index(const double eta, const double phi) const;
1607  void _tj_set_jetinfo ( TiledJet * const jet, const int _jets_index);
1608  void _bj_remove_from_tiles(TiledJet * const jet);
1609  void _initialise_tiles();
1610  void _print_tiles(TiledJet * briefjets ) const;
1611  void _add_neighbours_to_tile_union(const int tile_index,
1612  std::vector<int> & tile_union, int & n_near_tiles) const;
1613  void _add_untagged_neighbours_to_tile_union(const int tile_index,
1614  std::vector<int> & tile_union, int & n_near_tiles);
1615  struct EEBriefJet {
1616  double NN_dist; // obligatorily present
1617  double kt2; // obligatorily present == E^2 in general
1618  EEBriefJet * NN; // must be present too
1619  int _jets_index; // must also be present!
1620  double nx, ny, nz; // our internal storage for fast distance calcs
1621  };
1622  void _simple_N2_cluster_BriefJet();
1623  void _simple_N2_cluster_EEBriefJet();
1624 };
1625 template<class L> void ClusterSequence::_transfer_input_jets(
1626  const std::vector<L> & pseudojets) {
1627  _jets.reserve(pseudojets.size()*2);
1628  for (unsigned int i = 0; i < pseudojets.size(); i++) {
1629  _jets.push_back(pseudojets[i]);}
1630 }
1631 template<class L> ClusterSequence::ClusterSequence (
1632  const std::vector<L> & pseudojets,
1633  const JetDefinition & jet_def_in,
1634  const bool & writeout_combinations) :
1635  _jet_def(jet_def_in), _writeout_combinations(writeout_combinations),
1636  _structure_shared_ptr(new ClusterSequenceStructure(this))
1637 {
1638  _transfer_input_jets(pseudojets);
1639  _decant_options_partial();
1640  _initialise_and_run_no_decant();
1641 }
1642 inline const std::vector<PseudoJet> & ClusterSequence::jets () const {
1643  return _jets;
1644 }
1645 inline const std::vector<ClusterSequence::history_element> & ClusterSequence::history () const {
1646  return _history;
1647 }
1648 inline unsigned int ClusterSequence::n_particles() const {return _initial_n;}
1649 #ifndef __CINT__
1650 template<class L>
1651 std::vector<PseudoJet> JetDefinition::operator()(const std::vector<L> & particles) const {
1652  ClusterSequence * cs = new ClusterSequence(particles, *this);
1653  std::vector<PseudoJet> jets;
1654  if (is_spherical()) {
1655  jets = sorted_by_E(cs->inclusive_jets());
1656  } else {
1657  jets = sorted_by_pt(cs->inclusive_jets());
1658  }
1659  if (jets.size() != 0) {
1660  cs->delete_self_when_unused();
1661  } else {
1662  delete cs;
1663  }
1664  return jets;
1665 }
1666 #endif // __CINT__
1667 template <class J> inline void ClusterSequence::_bj_set_jetinfo(
1668  J * const jetA, const int _jets_index) const {
1669  jetA->eta = _jets[_jets_index].rap();
1670  jetA->phi = _jets[_jets_index].phi_02pi();
1671  jetA->kt2 = jet_scale_for_algorithm(_jets[_jets_index]);
1672  jetA->_jets_index = _jets_index;
1673  jetA->NN_dist = _R2;
1674  jetA->NN = NULL;
1675 }
1676 template <class J> inline double ClusterSequence::_bj_dist(
1677  const J * const jetA, const J * const jetB) const {
1678 #ifndef FJCORE_NEW_DELTA_PHI
1679  double dphi = std::abs(jetA->phi - jetB->phi);
1680  double deta = (jetA->eta - jetB->eta);
1681  if (dphi > pi) {dphi = twopi - dphi;}
1682 #else
1683  double dphi = pi-std::abs(pi-std::abs(jetA->phi - jetB->phi));
1684  double deta = (jetA->eta - jetB->eta);
1685 #endif
1686  return dphi*dphi + deta*deta;
1687 }
1688 template <class J> inline double ClusterSequence::_bj_diJ(const J * const jet) const {
1689  double kt2 = jet->kt2;
1690  if (jet->NN != NULL) {if (jet->NN->kt2 < kt2) {kt2 = jet->NN->kt2;}}
1691  return jet->NN_dist * kt2;
1692 }
1693 template <class J> inline void ClusterSequence::_bj_set_NN_nocross(
1694  J * const jet, J * const head, const J * const tail) const {
1695  double NN_dist = _R2;
1696  J * NN = NULL;
1697  if (head < jet) {
1698  for (J * jetB = head; jetB != jet; jetB++) {
1699  double dist = _bj_dist(jet,jetB);
1700  if (dist < NN_dist) {
1701  NN_dist = dist;
1702  NN = jetB;
1703  }
1704  }
1705  }
1706  if (tail > jet) {
1707  for (J * jetB = jet+1; jetB != tail; jetB++) {
1708  double dist = _bj_dist(jet,jetB);
1709  if (dist < NN_dist) {
1710  NN_dist = dist;
1711  NN = jetB;
1712  }
1713  }
1714  }
1715  jet->NN = NN;
1716  jet->NN_dist = NN_dist;
1717 }
1718 template <class J> inline void ClusterSequence::_bj_set_NN_crosscheck(J * const jet,
1719  J * const head, const J * const tail) const {
1720  double NN_dist = _R2;
1721  J * NN = NULL;
1722  for (J * jetB = head; jetB != tail; jetB++) {
1723  double dist = _bj_dist(jet,jetB);
1724  if (dist < NN_dist) {
1725  NN_dist = dist;
1726  NN = jetB;
1727  }
1728  if (dist < jetB->NN_dist) {
1729  jetB->NN_dist = dist;
1730  jetB->NN = jet;
1731  }
1732  }
1733  jet->NN = NN;
1734  jet->NN_dist = NN_dist;
1735 }
1736 FJCORE_END_NAMESPACE
1737 #endif // __FJCORE_CLUSTERSEQUENCE_HH__
1738 #ifndef __FJCORE_NNBASE_HH__
1739 #define __FJCORE_NNBASE_HH__
1740 FJCORE_BEGIN_NAMESPACE // defined in fastjet/internal/base.hh
1741 class _NoInfo {};
1742 template<class I> class NNInfo {
1743 public:
1744  NNInfo() : _info(NULL) {}
1745  NNInfo(I * info) : _info(info) {}
1746  template<class BJ> void init_jet(BJ * briefjet, const fjcore::PseudoJet & jet, int index) { briefjet->init(jet, index, _info);}
1747 private:
1748  I * _info;
1749 };
1750 template<> class NNInfo<_NoInfo> {
1751 public:
1752  NNInfo() {}
1753  NNInfo(_NoInfo * ) {}
1754  template<class BJ> void init_jet(BJ * briefjet, const fjcore::PseudoJet & jet, int index) { briefjet->init(jet, index);}
1755 };
1756 template<class I = _NoInfo> class NNBase : public NNInfo<I> {
1757 public:
1758  NNBase() {}
1759  NNBase(I * info) : NNInfo<I>(info) {}
1760  virtual void start(const std::vector<PseudoJet> & jets) = 0;
1761  virtual double dij_min(int & iA, int & iB) = 0;
1762  virtual void remove_jet(int iA) = 0;
1763  virtual void merge_jets(int iA, int iB, const PseudoJet & jet, int jet_index) = 0;
1764  virtual ~NNBase() {};
1765 };
1766 FJCORE_END_NAMESPACE // defined in fastjet/internal/base.hh
1767 #endif // __FJCORE_NNBASE_HH__
1768 #ifndef __FJCORE_NNH_HH__
1769 #define __FJCORE_NNH_HH__
1770 FJCORE_BEGIN_NAMESPACE // defined in fastjet/internal/base.hh
1771 template<class BJ, class I = _NoInfo> class NNH : public NNBase<I> {
1772 public:
1773  NNH(const std::vector<PseudoJet> & jets) : NNBase<I>() {start(jets);}
1774  NNH(const std::vector<PseudoJet> & jets, I * info) : NNBase<I>(info) {start(jets);}
1775  void start(const std::vector<PseudoJet> & jets);
1776  double dij_min(int & iA, int & iB);
1777  void remove_jet(int iA);
1778  void merge_jets(int iA, int iB, const PseudoJet & jet, int jet_index);
1779  ~NNH() {
1780  delete[] briefjets;
1781  }
1782 private:
1783  class NNBJ; // forward declaration
1784  void set_NN_crosscheck(NNBJ * jet, NNBJ * begin, NNBJ * end);
1785  void set_NN_nocross (NNBJ * jet, NNBJ * begin, NNBJ * end);
1786  NNBJ * briefjets;
1787  NNBJ * head, * tail;
1788  int n;
1789  std::vector<NNBJ *> where_is;
1790  class NNBJ : public BJ {
1791  public:
1792  void init(const PseudoJet & jet, int index_in) {
1793  BJ::init(jet);
1794  other_init(index_in);
1795  }
1796  void init(const PseudoJet & jet, int index_in, I * info) {
1797  BJ::init(jet, info);
1798  other_init(index_in);
1799  }
1800  void other_init(int index_in) {
1801  _index = index_in;
1802  NN_dist = BJ::beam_distance();
1803  NN = NULL;
1804  }
1805  int index() const {return _index;}
1806  double NN_dist;
1807  NNBJ * NN;
1808  private:
1809  int _index;
1810  };
1811 };
1812 template<class BJ, class I> void NNH<BJ,I>::start(const std::vector<PseudoJet> & jets) {
1813  n = jets.size();
1814  briefjets = new NNBJ[n];
1815  where_is.resize(2*n);
1816  NNBJ * jetA = briefjets;
1817  for (int i = 0; i< n; i++) {
1818  this->init_jet(jetA, jets[i], i);
1819  where_is[i] = jetA;
1820  jetA++; // move on to next entry of briefjets
1821  }
1822  tail = jetA; // a semaphore for the end of briefjets
1823  head = briefjets; // a nicer way of naming start
1824  for (jetA = head + 1; jetA != tail; jetA++) {
1825  set_NN_crosscheck(jetA, head, jetA);
1826  }
1827 }
1828 template<class BJ, class I> double NNH<BJ,I>::dij_min(int & iA, int & iB) {
1829  double diJ_min = briefjets[0].NN_dist;
1830  int diJ_min_jet = 0;
1831  for (int i = 1; i < n; i++) {
1832  if (briefjets[i].NN_dist < diJ_min) {
1833  diJ_min_jet = i;
1834  diJ_min = briefjets[i].NN_dist;
1835  }
1836  }
1837  NNBJ * jetA = & briefjets[diJ_min_jet];
1838  iA = jetA->index();
1839  iB = jetA->NN ? jetA->NN->index() : -1;
1840  return diJ_min;
1841 }
1842 template<class BJ, class I> void NNH<BJ,I>::remove_jet(int iA) {
1843  NNBJ * jetA = where_is[iA];
1844  tail--; n--;
1845  *jetA = *tail;
1846  where_is[jetA->index()] = jetA;
1847  for (NNBJ * jetI = head; jetI != tail; jetI++) {
1848  if (jetI->NN == jetA) set_NN_nocross(jetI, head, tail);
1849  if (jetI->NN == tail) {jetI->NN = jetA;}
1850  }
1851 }
1852 template<class BJ, class I> void NNH<BJ,I>::merge_jets(int iA, int iB,
1853  const PseudoJet & jet, int index) {
1854  NNBJ * jetA = where_is[iA];
1855  NNBJ * jetB = where_is[iB];
1856  if (jetA < jetB) std::swap(jetA,jetB);
1857  this->init_jet(jetB, jet, index);
1858  if (index >= int(where_is.size())) where_is.resize(2*index);
1859  where_is[jetB->index()] = jetB;
1860  tail--; n--;
1861  *jetA = *tail;
1862  where_is[jetA->index()] = jetA;
1863  for (NNBJ * jetI = head; jetI != tail; jetI++) {
1864  if (jetI->NN == jetA || jetI->NN == jetB) {
1865  set_NN_nocross(jetI, head, tail);
1866  }
1867  double dist = jetI->distance(jetB);
1868  if (dist < jetI->NN_dist) {
1869  if (jetI != jetB) {
1870  jetI->NN_dist = dist;
1871  jetI->NN = jetB;
1872  }
1873  }
1874  if (dist < jetB->NN_dist) {
1875  if (jetI != jetB) {
1876  jetB->NN_dist = dist;
1877  jetB->NN = jetI;
1878  }
1879  }
1880  if (jetI->NN == tail) {jetI->NN = jetA;}
1881  }
1882 }
1883 template <class BJ, class I> void NNH<BJ,I>::set_NN_crosscheck(NNBJ * jet,
1884  NNBJ * begin, NNBJ * end) {
1885  double NN_dist = jet->beam_distance();
1886  NNBJ * NN = NULL;
1887  for (NNBJ * jetB = begin; jetB != end; jetB++) {
1888  double dist = jet->distance(jetB);
1889  if (dist < NN_dist) {
1890  NN_dist = dist;
1891  NN = jetB;
1892  }
1893  if (dist < jetB->NN_dist) {
1894  jetB->NN_dist = dist;
1895  jetB->NN = jet;
1896  }
1897  }
1898  jet->NN = NN;
1899  jet->NN_dist = NN_dist;
1900 }
1901 template <class BJ, class I> void NNH<BJ,I>::set_NN_nocross(
1902  NNBJ * jet, NNBJ * begin, NNBJ * end) {
1903  double NN_dist = jet->beam_distance();
1904  NNBJ * NN = NULL;
1905  if (begin < jet) {
1906  for (NNBJ * jetB = begin; jetB != jet; jetB++) {
1907  double dist = jet->distance(jetB);
1908  if (dist < NN_dist) {
1909  NN_dist = dist;
1910  NN = jetB;
1911  }
1912  }
1913  }
1914  if (end > jet) {
1915  for (NNBJ * jetB = jet+1; jetB != end; jetB++) {
1916  double dist = jet->distance (jetB);
1917  if (dist < NN_dist) {
1918  NN_dist = dist;
1919  NN = jetB;
1920  }
1921  }
1922  }
1923  jet->NN = NN;
1924  jet->NN_dist = NN_dist;
1925 }
1926 FJCORE_END_NAMESPACE // defined in fastjet/internal/base.hh
1927 #endif // __FJCORE_NNH_HH__
1928 #endif
Error
Definition: FJcore.h:367
integral_type
Definition: FJcore.h:100
NNH
Definition: FJcore.h:1604
integral_type::value
static const T value
the value (only member carrying info)
Definition: FJcore.h:101
integral_type::value_type
T value_type
a typedef for the type T
Definition: FJcore.h:200
CompositeJetStructure::_pieces
std::vector< PseudoJet > _pieces
the pieces building the jet
Definition: FJcore.h:1112
CompositeJetStructure::_area_4vector_ptr
PseudoJet * _area_4vector_ptr
pointer to the 4-vector jet area
Definition: FJcore.h:1113
Tile
Definition: FJcore.cc:837
__inheritance_helper
Definition: FJcore.h:112
_NoInfo
Definition: FJcore.h:1589
integral_type::type
integral_type< T, _t > type
a typedef for the whole structure
Definition: FJcore.h:201
Generated by   doxygen 1.8.5