Analysis Software
Documentation for sPHENIX simulation software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
fjcore.hh
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file fjcore.hh
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 #define FJCORE_BEGIN_NAMESPACE namespace fjcore {
172 #define FJCORE_END_NAMESPACE }
173 #ifdef FJCORE_HAVE_OVERRIDE
174 # define FJCORE_OVERRIDE override
175 #else
176 # define FJCORE_OVERRIDE
177 #endif
178 #endif // __FJCORE_FASTJET_BASE_HH__
179 #ifndef __FJCORE_NUMCONSTS__
180 #define __FJCORE_NUMCONSTS__
181 FJCORE_BEGIN_NAMESPACE // defined in fastjet/internal/base.hh
182 const double pi = 3.141592653589793238462643383279502884197;
183 const double twopi = 6.283185307179586476925286766559005768394;
184 const double pisq = 9.869604401089358618834490999876151135314;
185 const double zeta2 = 1.644934066848226436472415166646025189219;
186 const double zeta3 = 1.202056903159594285399738161511449990765;
187 const double eulergamma = 0.577215664901532860606512090082402431042;
188 const double ln2 = 0.693147180559945309417232121458176568076;
190 #endif // __FJCORE_NUMCONSTS__
191 #ifndef __FJCORE_INTERNAL_IS_BASE_HH__
192 #define __FJCORE_INTERNAL_IS_BASE_HH__
194 template<typename T, T _t>
196  static const T value = _t;
197  typedef T value_type;
199 };
200 template<typename T, T _t>
204 typedef char (&__yes_type)[1]; //< the yes type
205 typedef char (&__no_type) [2]; //< the no type
206 template<typename B, typename D>
208 #if !((_MSC_VER !=0 ) && (_MSC_VER == 1310)) // MSVC 7.1
209  template <typename T>
210  static __yes_type check_sig(D const volatile *, T);
211 #else
212  static __yes_type check_sig(D const volatile *, long);
213 #endif
214  static __no_type check_sig(B const volatile *, int);
215 };
216 template<typename B, typename D>
218 #if ((_MSC_FULL_VER != 0) && (_MSC_FULL_VER >= 140050000))
219 #pragma warning(push)
220 #pragma warning(disable:6334)
221 #endif
222  struct Host{
223 #if !((_MSC_VER !=0 ) && (_MSC_VER == 1310))
224  operator B const volatile *() const;
225 #else
226  operator B const volatile * const&() const;
227 #endif
228  operator D const volatile *();
229  };
230  static const bool value = ((sizeof(B)!=0) &&
231  (sizeof(D)!=0) &&
232  (sizeof(__inheritance_helper<B,D>::check_sig(Host(), 0)) == sizeof(__yes_type)));
233 #if ((_MSC_FULL_VER != 0) && (_MSC_FULL_VER >= 140050000))
234 #pragma warning(pop)
235 #endif
236 };
237 template<class B, class D>
239  return IsBaseAndDerived<B,D>::value ? (B*)(d) : 0;
240 }
242 #endif // __IS_BASE_OF_HH__
243 #ifndef __FJCORE_FJCORE_DEPRECATED_HH__
244 #define __FJCORE_FJCORE_DEPRECATED_HH__
245 #if defined(FJCORE_HAVE_CXX14_DEPRECATED)
246 # define FJCORE_DEPRECATED [[deprecated]]
247 # define FJCORE_DEPRECATED_MSG(message) [[deprecated(message)]]
248 #elif defined(FJCORE_HAVE_GNUCXX_DEPRECATED)
249 # define FJCORE_DEPRECATED __attribute__((__deprecated__))
250 # define FJCORE_DEPRECATED_MSG(message) __attribute__((__deprecated__))
251 #else
252 # define FJCORE_DEPRECATED
253 # define FJCORE_DEPRECATED_MSG(message)
254 #endif
255 #endif // __FJCORE_FJCORE_DEPRECATED_HH__
256 #ifndef __FJCORE_SHARED_PTR_HH__
257 #define __FJCORE_SHARED_PTR_HH__
258 #include <cstdlib> // for NULL!!!
259 #ifdef __FJCORE_USETR1SHAREDPTR
260 #include <tr1/memory>
261 #endif // __FJCORE_USETR1SHAREDPTR
262 FJCORE_BEGIN_NAMESPACE // defined in fastjet/internal/base.hh
263 #ifdef __FJCORE_USETR1SHAREDPTR
264 template<class T>
265 class SharedPtr : public std::tr1::shared_ptr<T> {
266 public:
267  SharedPtr() : std::tr1::shared_ptr<T>() {}
268  SharedPtr(T * t) : std::tr1::shared_ptr<T>(t) {}
269  SharedPtr(const SharedPtr<T> & t) : std::tr1::shared_ptr<T>(t) {}
270  #ifdef FJCORE_HAVE_EXPLICIT_FOR_OPERATORS
271  explicit
272  #endif
273  inline operator bool() const {return (this->get()!=NULL);}
274  T* operator ()() const{
275  return this->get(); // automatically returns NULL when out-of-scope
276  }
277 };
278 #else // __FJCORE_USETR1SHAREDPTR
279 template<class T>
280 class SharedPtr{
281 public:
282  class __SharedCountingPtr;
283  SharedPtr() : _ptr(NULL){}
284  template<class Y> explicit SharedPtr(Y* ptr){
285  _ptr = new __SharedCountingPtr(ptr);
286  }
287  SharedPtr(SharedPtr const & share) : _ptr(share._get_container()){
288  if (_ptr!=NULL) ++(*_ptr);
289  }
291  if (_ptr==NULL) return;
292  _decrease_count();
293  }
294  void reset(){
295  SharedPtr().swap(*this);
296  }
297  template<class Y> void reset(Y * ptr){
298  SharedPtr(ptr).swap(*this);
299  }
300  template<class Y> void reset(SharedPtr<Y> const & share){
301  if (_ptr!=NULL){
302  if (_ptr == share._get_container()) return;
303  _decrease_count();
304  }
305  _ptr = share._get_container(); // Note: automatically set it to NULL if share is empty
306  if (_ptr!=NULL) ++(*_ptr);
307  }
308  SharedPtr& operator=(SharedPtr const & share){
309  reset(share);
310  return *this;
311  }
312  template<class Y> SharedPtr& operator=(SharedPtr<Y> const & share){
313  reset(share);
314  return *this;
315  }
316  FJCORE_DEPRECATED_MSG("Use SharedPtr<T>::get() instead")
317  T* operator ()() const{
318  if (_ptr==NULL) return NULL;
319  return _ptr->get(); // automatically returns NULL when out-of-scope
320  }
321  inline T& operator*() const{
322  return *(_ptr->get());
323  }
324  inline T* operator->() const{
325  if (_ptr==NULL) return NULL;
326  return _ptr->get();
327  }
328  inline T* get() const{
329  if (_ptr==NULL) return NULL;
330  return _ptr->get();
331  }
332  inline bool unique() const{
333  return (use_count()==1);
334  }
335  inline long use_count() const{
336  if (_ptr==NULL) return 0;
337  return _ptr->use_count(); // automatically returns NULL when out-of-scope
338  }
339  #ifdef FJCORE_HAVE_EXPLICIT_FOR_OPERATORS
340  explicit
341  #endif
342  inline operator bool() const{
343  return (get()!=NULL);
344  }
345  inline void swap(SharedPtr & share){
346  __SharedCountingPtr* share_container = share._ptr;
347  share._ptr = _ptr;
348  _ptr = share_container;
349  }
350  void set_count(const long & count){
351  if (_ptr==NULL) return;
352  _ptr->set_count(count);
353  }
355  public:
357  template<class Y> explicit __SharedCountingPtr(Y* ptr) : _ptr(ptr), _count(1){}
359  if (_ptr!=NULL){ delete _ptr;}
360  }
361  inline T* get() const {return _ptr;}
362  inline long use_count() const {return _count;}
363  inline long operator++(){return ++_count;}
364  inline long operator--(){return --_count;}
365  inline long operator++(int){return _count++;}
366  inline long operator--(int){return _count--;}
367  void set_count(const long & count){
368  _count = count;
369  }
370  private:
371  T *_ptr;
372  long _count;
373  };
374 private:
375  inline __SharedCountingPtr* _get_container() const{
376  return _ptr;
377  }
379  (*_ptr)--;
380  if (_ptr->use_count()==0)
381  delete _ptr; // that automatically deletes the object itself
382  }
383  __SharedCountingPtr *_ptr;
384 };
385 template<class T,class U>
386 inline bool operator==(SharedPtr<T> const & t, SharedPtr<U> const & u){
387  return t.get() == u.get();
388 }
389 template<class T,class U>
390 inline bool operator!=(SharedPtr<T> const & t, SharedPtr<U> const & u){
391  return t.get() != u.get();
392 }
393 template<class T,class U>
394 inline bool operator<(SharedPtr<T> const & t, SharedPtr<U> const & u){
395  return t.get() < u.get();
396 }
397 template<class T>
398 inline void swap(SharedPtr<T> & a, SharedPtr<T> & b){
399  return a.swap(b);
400 }
401 template<class T>
402 inline T* get_pointer(SharedPtr<T> const & t){
403  return t.get();
404 }
405 #endif // __FJCORE_USETR1SHAREDPTR
406 FJCORE_END_NAMESPACE // defined in fastjet/internal/base.hh
407 #endif // __FJCORE_SHARED_PTR_HH__
408 #ifndef __FJCORE_LIMITEDWARNING_HH__
409 #define __FJCORE_LIMITEDWARNING_HH__
410 #include <iostream>
411 #include <string>
412 #include <list>
415 public:
417  LimitedWarning(int max_warn_in) : _max_warn(max_warn_in), _n_warn_so_far(0), _this_warning_summary(0) {}
418  void warn(const char * warning) {warn(warning, _default_ostr);}
419  void warn(const std::string & warning) {warn(warning.c_str(), _default_ostr);}
420  void warn(const char * warning, std::ostream * ostr);
421  void warn(const std::string & warning, std::ostream * ostr) {warn(warning.c_str(), ostr);}
422  static void set_default_stream(std::ostream * ostr) {
423  _default_ostr = ostr;
424  }
425  static void set_default_max_warn(int max_warn) {
427  }
428  int max_warn() const {return _max_warn;}
429  int n_warn_so_far() const {return _n_warn_so_far;}
430  static std::string summary();
431 private:
433  static int _max_warn_default;
434  static std::ostream * _default_ostr;
435  typedef std::pair<std::string, unsigned int> Summary;
436  static std::list< Summary > _global_warnings_summary;
438 };
440 #endif // __FJCORE_LIMITEDWARNING_HH__
441 #ifndef __FJCORE_ERROR_HH__
442 #define __FJCORE_ERROR_HH__
443 #include<iostream>
444 #include<string>
445 #if (!defined(FJCORE_HAVE_EXECINFO_H)) || defined(__FJCORE__)
446 #endif
447 FJCORE_BEGIN_NAMESPACE // defined in fastjet/internal/base.hh
448 class Error {
449 public:
450  Error() {}
451  Error(const std::string & message);
452  virtual ~Error() {}
453  std::string message() const {return _message;}
454  static void set_print_errors(bool print_errors) {_print_errors = print_errors;}
455  static void set_print_backtrace(bool enabled);
456  static void set_default_stream(std::ostream * ostr) {
457  _default_ostr = ostr;
458  }
459 private:
461  static bool _print_errors;
462  static bool _print_backtrace;
463  static std::ostream * _default_ostr;
464 #if (!defined(FJCORE_HAVE_EXECINFO_H)) || defined(__FJCORE__)
466 #endif
467 };
468 class InternalError : public Error{
469 public:
470  InternalError(const std::string & message_in) : Error(std::string("*** CRITICAL INTERNAL FASTJET ERROR *** CONTACT THE AUTHORS *** ") + message_in){ }
471 };
473 #endif // __FJCORE_ERROR_HH__
474 #ifndef __FJCORE_PSEUDOJET_STRUCTURE_BASE_HH__
475 #define __FJCORE_PSEUDOJET_STRUCTURE_BASE_HH__
476 #include <vector>
477 #include <string>
478 FJCORE_BEGIN_NAMESPACE // defined in fastjet/internal/base.hh
479 class PseudoJet;
480 class ClusterSequence;
482 public:
485  virtual std::string description() const{ return "PseudoJet with an unknown structure"; }
486  virtual bool has_associated_cluster_sequence() const { return false;}
487  virtual const ClusterSequence* associated_cluster_sequence() const;
488  virtual bool has_valid_cluster_sequence() const {return false;}
489  virtual const ClusterSequence * validated_cs() const;
490  virtual bool has_partner(const PseudoJet &reference, PseudoJet &partner) const;
491  virtual bool has_child(const PseudoJet &reference, PseudoJet &child) const;
492  virtual bool has_parents(const PseudoJet &reference, PseudoJet &parent1, PseudoJet &parent2) const;
493  virtual bool object_in_jet(const PseudoJet &reference, const PseudoJet &jet) const;
494  virtual bool has_constituents() const {return false;}
495  virtual std::vector<PseudoJet> constituents(const PseudoJet &reference) const;
496  virtual bool has_exclusive_subjets() const {return false;}
497  virtual std::vector<PseudoJet> exclusive_subjets(const PseudoJet &reference, const double & dcut) const;
498  virtual int n_exclusive_subjets(const PseudoJet &reference, const double & dcut) const;
499  virtual std::vector<PseudoJet> exclusive_subjets_up_to (const PseudoJet &reference, int nsub) const;
500  virtual double exclusive_subdmerge(const PseudoJet &reference, int nsub) const;
501  virtual double exclusive_subdmerge_max(const PseudoJet &reference, int nsub) const;
502  virtual bool has_pieces(const PseudoJet & /* reference */) const {
503  return false;}
504  virtual std::vector<PseudoJet> pieces(const PseudoJet & /* reference */
505  ) const;
506 };
508 #endif // __FJCORE_PSEUDOJET_STRUCTURE_BASE_HH__
509 #ifndef __FJCORE_PSEUDOJET_HH__
510 #define __FJCORE_PSEUDOJET_HH__
511 #include<valarray>
512 #include<vector>
513 #include<cassert>
514 #include<cmath>
515 #include<iostream>
516 FJCORE_BEGIN_NAMESPACE // defined in fastjet/internal/base.hh
517 const double MaxRap = 1e5;
518 const double pseudojet_invalid_phi = -100.0;
519 const double pseudojet_invalid_rap = -1e200;
520 class PseudoJet {
521  public:
522  PseudoJet() : _px(0), _py(0), _pz(0), _E(0) {_finish_init(); _reset_indices();}
523  PseudoJet(const double px, const double py, const double pz, const double E);
524  template <class L> PseudoJet(const L & some_four_vector);
525  PseudoJet(bool /* dummy */) {}
526  virtual ~PseudoJet(){};
527  inline double E() const {return _E;}
528  inline double e() const {return _E;} // like CLHEP
529  inline double px() const {return _px;}
530  inline double py() const {return _py;}
531  inline double pz() const {return _pz;}
532  inline double phi() const {return phi_02pi();}
533  inline double phi_std() const {
535  return _phi > pi ? _phi-twopi : _phi;}
536  inline double phi_02pi() const {
538  return _phi;
539  }
540  inline double rap() const {
542  return _rap;
543  }
544  inline double rapidity() const {return rap();} // like CLHEP
545  double pseudorapidity() const;
546  double eta() const {return pseudorapidity();}
547  inline double pt2() const {return _kt2;}
548  inline double pt() const {return sqrt(_kt2);}
549  inline double perp2() const {return _kt2;} // like CLHEP
550  inline double perp() const {return sqrt(_kt2);} // like CLHEP
551  inline double kt2() const {return _kt2;} // for bkwds compatibility
552  inline double m2() const {return (_E+_pz)*(_E-_pz)-_kt2;}
553  inline double m() const;
554  inline double mperp2() const {return (_E+_pz)*(_E-_pz);}
555  inline double mperp() const {return sqrt(std::abs(mperp2()));}
556  inline double mt2() const {return (_E+_pz)*(_E-_pz);}
557  inline double mt() const {return sqrt(std::abs(mperp2()));}
558  inline double modp2() const {return _kt2+_pz*_pz;}
559  inline double modp() const {return sqrt(_kt2+_pz*_pz);}
560  inline double Et() const {return (_kt2==0) ? 0.0 : _E/sqrt(1.0+_pz*_pz/_kt2);}
561  inline double Et2() const {return (_kt2==0) ? 0.0 : _E*_E/(1.0+_pz*_pz/_kt2);}
562  double operator () (int i) const ;
563  inline double operator [] (int i) const { return (*this)(i); }; // this too
564  double kt_distance(const PseudoJet & other) const;
565  double plain_distance(const PseudoJet & other) const;
566  inline double squared_distance(const PseudoJet & other) const {
567  return plain_distance(other);}
568  inline double delta_R(const PseudoJet & other) const {
569  return sqrt(squared_distance(other));
570  }
571  double delta_phi_to(const PseudoJet & other) const;
572  inline double beam_distance() const {return _kt2;}
573  std::valarray<double> four_mom() const;
574  enum { X=0, Y=1, Z=2, T=3, NUM_COORDINATES=4, SIZE=NUM_COORDINATES };
575  PseudoJet & boost(const PseudoJet & prest);
576  PseudoJet & unboost(const PseudoJet & prest);
577  void operator*=(double);
578  void operator/=(double);
579  void operator+=(const PseudoJet &);
580  void operator-=(const PseudoJet &);
581  inline void reset(double px, double py, double pz, double E);
582  inline void reset(const PseudoJet & psjet) {
583  (*this) = psjet;
584  }
585  template <class L> inline void reset(const L & some_four_vector) {
586  const PseudoJet * pj = fjcore::cast_if_derived<const PseudoJet>(&some_four_vector);
587  if (pj){
588  (*this) = *pj;
589  } else {
590  reset(some_four_vector[0], some_four_vector[1],
591  some_four_vector[2], some_four_vector[3]);
592  }
593  }
594  inline void reset_PtYPhiM(double pt_in, double y_in, double phi_in, double m_in=0.0) {
595  reset_momentum_PtYPhiM(pt_in, y_in, phi_in, m_in);
596  _reset_indices();
597  }
598  inline void reset_momentum(double px, double py, double pz, double E);
599  inline void reset_momentum(const PseudoJet & pj);
600  void reset_momentum_PtYPhiM(double pt, double y, double phi, double m=0.0);
601  template <class L> inline void reset_momentum(const L & some_four_vector) {
602  reset_momentum(some_four_vector[0], some_four_vector[1],
603  some_four_vector[2], some_four_vector[3]);
604  }
605  void set_cached_rap_phi(double rap, double phi);
606  inline int user_index() const {return _user_index;}
607  inline void set_user_index(const int index) {_user_index = index;}
609  public:
611  virtual ~UserInfoBase(){};
612  };
613  class InexistentUserInfo : public Error {
614  public:
616  };
617  void set_user_info(UserInfoBase * user_info_in) {
618  _user_info.reset(user_info_in);
619  }
620  template<class L>
621  const L & user_info() const{
622  if (_user_info.get() == 0) throw InexistentUserInfo();
623  return dynamic_cast<const L &>(* _user_info.get());
624  }
625  bool has_user_info() const{
626  return _user_info.get();
627  }
628  template<class L>
629  bool has_user_info() const{
630  return _user_info.get() && dynamic_cast<const L *>(_user_info.get());
631  }
632  const UserInfoBase * user_info_ptr() const{
633  return _user_info.get();
634  }
636  return _user_info;
637  }
639  return _user_info;
640  }
641  std::string description() const;
642  bool has_associated_cluster_sequence() const;
644  bool has_valid_cluster_sequence() const;
645  bool has_valid_cs() const {return has_valid_cluster_sequence();}
649  return validated_cs();
650  }
651  const ClusterSequence * validated_cs() const;
653  bool has_structure() const;
654  const PseudoJetStructureBase* structure_ptr() const;
658  template<typename StructureType>
659  const StructureType & structure() const;
660  template<typename TransformerType>
661  bool has_structure_of() const;
662  template<typename TransformerType>
663  const typename TransformerType::StructureType & structure_of() const;
664  virtual bool has_partner(PseudoJet &partner) const;
665  virtual bool has_child(PseudoJet &child) const;
666  virtual bool has_parents(PseudoJet &parent1, PseudoJet &parent2) const;
667  virtual bool contains(const PseudoJet &constituent) const;
668  virtual bool is_inside(const PseudoJet &jet) const;
669  virtual bool has_constituents() const;
670  virtual std::vector<PseudoJet> constituents() const;
671  virtual bool has_exclusive_subjets() const;
672  std::vector<PseudoJet> exclusive_subjets (const double dcut) const;
673  int n_exclusive_subjets(const double dcut) const;
674  std::vector<PseudoJet> exclusive_subjets (int nsub) const;
675  std::vector<PseudoJet> exclusive_subjets_up_to (int nsub) const;
676  double exclusive_subdmerge(int nsub) const;
677  double exclusive_subdmerge_max(int nsub) const;
678  virtual bool has_pieces() const;
679  virtual std::vector<PseudoJet> pieces() const;
680  inline int cluster_hist_index() const {return _cluster_hist_index;}
682  inline int cluster_sequence_history_index() const {
683  return cluster_hist_index();}
684  inline void set_cluster_sequence_history_index(const int index) {
685  set_cluster_hist_index(index);}
686  protected:
689  private:
690  double _px,_py,_pz,_E;
691  mutable double _phi, _rap;
692  double _kt2;
694  void _finish_init();
695  void _reset_indices();
696  inline void _ensure_valid_rap_phi() const {
697  if (_phi == pseudojet_invalid_phi) _set_rap_phi();
698  }
699  void _set_rap_phi() const;
700  friend PseudoJet operator*(double, const PseudoJet &);
701 };
702 PseudoJet operator+(const PseudoJet &, const PseudoJet &);
703 PseudoJet operator-(const PseudoJet &, const PseudoJet &);
704 PseudoJet operator*(double, const PseudoJet &);
705 PseudoJet operator*(const PseudoJet &, double);
706 PseudoJet operator/(const PseudoJet &, double);
707 bool operator==(const PseudoJet &, const PseudoJet &);
708 inline bool operator!=(const PseudoJet & a, const PseudoJet & b) {return !(a==b);}
709 bool operator==(const PseudoJet & jet, const double val);
710 inline bool operator==(const double val, const PseudoJet & jet) {return jet == val;}
711 inline bool operator!=(const PseudoJet & a, const double val) {return !(a==val);}
712 inline bool operator!=( const double val, const PseudoJet & a) {return !(a==val);}
713 inline double dot_product(const PseudoJet & a, const PseudoJet & b) {
714  return a.E()*b.E() - a.px()*b.px() - a.py()*b.py() - a.pz()*b.pz();
715 }
716 bool have_same_momentum(const PseudoJet &, const PseudoJet &);
717 PseudoJet PtYPhiM(double pt, double y, double phi, double m = 0.0);
718 std::vector<PseudoJet> sorted_by_pt(const std::vector<PseudoJet> & jets);
719 std::vector<PseudoJet> sorted_by_rapidity(const std::vector<PseudoJet> & jets);
720 std::vector<PseudoJet> sorted_by_E(const std::vector<PseudoJet> & jets);
721 std::vector<PseudoJet> sorted_by_pz(const std::vector<PseudoJet> & jets);
722 void sort_indices(std::vector<int> & indices,
723  const std::vector<double> & values);
724 template<class T> std::vector<T> objects_sorted_by_values(const std::vector<T> & objects,
725  const std::vector<double> & values) {
726  if (objects.size() != values.size()){
727  throw Error("fjcore::objects_sorted_by_values(...): the size of the 'objects' vector must match the size of the 'values' vector");
728  }
729  std::vector<int> indices(values.size());
730  for (size_t i = 0; i < indices.size(); i++) {indices[i] = i;}
731  sort_indices(indices, values);
732  std::vector<T> objects_sorted(objects.size());
733  for (size_t i = 0; i < indices.size(); i++) {
734  objects_sorted[i] = objects[indices[i]];
735  }
736  return objects_sorted;
737 }
739 public:
740  inline IndexedSortHelper (const std::vector<double> * reference_values) {
741  _ref_values = reference_values;
742  };
743  inline int operator() (const int i1, const int i2) const {
744  return (*_ref_values)[i1] < (*_ref_values)[i2];
745  };
746 private:
747  const std::vector<double> * _ref_values;
748 };
749 template <class L> inline PseudoJet::PseudoJet(const L & some_four_vector) {
750  reset(some_four_vector);
751 }
754  set_user_index(-1);
755  _structure.reset();
756  _user_info.reset();
757 }
758 inline double PseudoJet::m() const {
759  double mm = m2();
760  return mm < 0.0 ? -std::sqrt(-mm) : std::sqrt(mm);
761 }
762 inline void PseudoJet::reset(double px_in, double py_in, double pz_in, double E_in) {
763  _px = px_in;
764  _py = py_in;
765  _pz = pz_in;
766  _E = E_in;
767  _finish_init();
768  _reset_indices();
769 }
770 inline void PseudoJet::reset_momentum(double px_in, double py_in, double pz_in, double E_in) {
771  _px = px_in;
772  _py = py_in;
773  _pz = pz_in;
774  _E = E_in;
775  _finish_init();
776 }
777 inline void PseudoJet::reset_momentum(const PseudoJet & pj) {
778  _px = pj._px ;
779  _py = pj._py ;
780  _pz = pj._pz ;
781  _E = pj._E ;
782  _phi = pj._phi;
783  _rap = pj._rap;
784  _kt2 = pj._kt2;
785 }
786 template<typename StructureType>
787 const StructureType & PseudoJet::structure() const{
788  return dynamic_cast<const StructureType &>(* validated_structure_ptr());
789 }
790 template<typename TransformerType>
792  if (!_structure) return false;
793  return dynamic_cast<const typename TransformerType::StructureType *>(_structure.get()) != 0;
794 }
795 template<typename TransformerType>
796 const typename TransformerType::StructureType & PseudoJet::structure_of() const{
797  if (!_structure)
798  throw Error("Trying to access the structure of a PseudoJet without an associated structure");
799  return dynamic_cast<const typename TransformerType::StructureType &>(*_structure);
800 }
801 PseudoJet join(const std::vector<PseudoJet> & pieces);
802 PseudoJet join(const PseudoJet & j1);
803 PseudoJet join(const PseudoJet & j1, const PseudoJet & j2);
804 PseudoJet join(const PseudoJet & j1, const PseudoJet & j2, const PseudoJet & j3);
805 PseudoJet join(const PseudoJet & j1, const PseudoJet & j2, const PseudoJet & j3, const PseudoJet & j4);
807 #endif // __FJCORE_PSEUDOJET_HH__
808 #ifndef __FJCORE_FUNCTION_OF_PSEUDOJET_HH__
809 #define __FJCORE_FUNCTION_OF_PSEUDOJET_HH__
811 template<typename TOut>
813 public:
816  virtual std::string description() const{ return "";}
817  virtual TOut result(const PseudoJet &pj) const = 0;
818  TOut operator()(const PseudoJet &pj) const { return result(pj);}
819  std::vector<TOut> operator()(const std::vector<PseudoJet> &pjs) const {
820  std::vector<TOut> res(pjs.size());
821  for (unsigned int i=0; i<pjs.size(); i++)
822  res[i] = result(pjs[i]);
823  return res;
824  }
825 };
827 #endif // __FJCORE_FUNCTION_OF_PSEUDOJET_HH__
828 #ifndef __FJCORE_SELECTOR_HH__
829 #define __FJCORE_SELECTOR_HH__
830 #include <limits>
831 #include <cmath>
832 FJCORE_BEGIN_NAMESPACE // defined in fastjet/internal/base.hh
833 class Selector;
835 public:
836  virtual ~SelectorWorker() {}
837  virtual bool pass(const PseudoJet & jet) const = 0;
838  virtual void terminator(std::vector<const PseudoJet *> & jets) const {
839  for (unsigned i = 0; i < jets.size(); i++) {
840  if (jets[i] && !pass(*jets[i])) jets[i] = NULL;
841  }
842  }
843  virtual bool applies_jet_by_jet() const {return true;}
844  virtual std::string description() const {return "missing description";}
845  virtual bool takes_reference() const { return false;}
846  virtual void set_reference(const PseudoJet & /*reference*/){
847  throw Error("set_reference(...) cannot be used for a selector worker that does not take a reference");
848  }
849  virtual SelectorWorker* copy(){
850  throw Error("this SelectorWorker has nothing to copy");
851  }
852  virtual void get_rapidity_extent(double & rapmin, double & rapmax) const {
853  rapmax = std::numeric_limits<double>::infinity();
854  rapmin = -rapmax;
855  }
856  virtual bool is_geometric() const { return false;}
857  virtual bool has_finite_area() const;
858  virtual bool has_known_area() const { return false;}
859  virtual double known_area() const{
860  throw Error("this selector has no computable area");
861  }
862 };
863 class Selector{
864 public:
865  Selector() {}
866  Selector(SelectorWorker * worker_in) {_worker.reset(worker_in);}
867  virtual ~Selector(){}
868  bool pass(const PseudoJet & jet) const {
870  throw Error("Cannot apply this selector to an individual jet");
871  }
872  return _worker->pass(jet);
873  }
874  bool operator()(const PseudoJet & jet) const {
875  return pass(jet);
876  }
877  unsigned int count(const std::vector<PseudoJet> & jets) const;
878  PseudoJet sum(const std::vector<PseudoJet> & jets) const;
879  double scalar_pt_sum(const std::vector<PseudoJet> & jets) const;
880  void sift(const std::vector<PseudoJet> & jets,
881  std::vector<PseudoJet> & jets_that_pass,
882  std::vector<PseudoJet> & jets_that_fail) const;
883  bool applies_jet_by_jet() const {
885  }
886  std::vector<PseudoJet> operator()(const std::vector<PseudoJet> & jets) const;
887  virtual void nullify_non_selected(std::vector<const PseudoJet *> & jets) const {
888  validated_worker()->terminator(jets);
889  }
890  void get_rapidity_extent(double &rapmin, double &rapmax) const {
891  return validated_worker()->get_rapidity_extent(rapmin, rapmax);
892  }
894  return validated_worker()->description();
895  }
896  bool is_geometric() const{
897  return validated_worker()->is_geometric();
898  }
899  bool has_finite_area() const{
900  return validated_worker()->has_finite_area();
901  }
902  const SharedPtr<SelectorWorker> & worker() const {return _worker;}
904  const SelectorWorker* worker_ptr = _worker.get();
905  if (worker_ptr == 0) throw InvalidWorker();
906  return worker_ptr;
907  }
908  bool takes_reference() const {
909  return validated_worker()->takes_reference();
910  }
911  const Selector & set_reference(const PseudoJet &reference){
912  if (! validated_worker()->takes_reference()){
913  return *this;
914  }
916  _worker->set_reference(reference);
917  return *this;
918  }
919  class InvalidWorker : public Error {
920  public:
921  InvalidWorker() : Error("Attempt to use Selector with no valid underlying worker") {}
922  };
923  class InvalidArea : public Error {
924  public:
925  InvalidArea() : Error("Attempt to obtain area from Selector for which this is not meaningful") {}
926  };
927  Selector & operator &=(const Selector & b);
928  Selector & operator |=(const Selector & b);
929 protected:
931  if (_worker.unique()) return;
933  }
934 private:
936 };
938 Selector operator!(const Selector & s);
939 Selector operator ||(const Selector & s1, const Selector & s2);
940 Selector operator&&(const Selector & s1, const Selector & s2);
941 Selector operator*(const Selector & s1, const Selector & s2);
942 Selector SelectorPtMin(double ptmin);
943 Selector SelectorPtMax(double ptmax);
944 Selector SelectorPtRange(double ptmin, double ptmax);
945 Selector SelectorEtMin(double Etmin);
946 Selector SelectorEtMax(double Etmax);
947 Selector SelectorEtRange(double Etmin, double Etmax);
948 Selector SelectorEMin(double Emin);
949 Selector SelectorEMax(double Emax);
950 Selector SelectorERange(double Emin, double Emax);
951 Selector SelectorMassMin(double Mmin);
952 Selector SelectorMassMax(double Mmax);
953 Selector SelectorMassRange(double Mmin, double Mmax);
954 Selector SelectorRapMin(double rapmin);
955 Selector SelectorRapMax(double rapmax);
956 Selector SelectorRapRange(double rapmin, double rapmax);
957 Selector SelectorAbsRapMin(double absrapmin);
958 Selector SelectorAbsRapMax(double absrapmax);
959 Selector SelectorAbsRapRange(double absrapmin, double absrapmax);
962 Selector SelectorEtaRange(double etamin, double etamax);
963 Selector SelectorAbsEtaMin(double absetamin);
964 Selector SelectorAbsEtaMax(double absetamax);
965 Selector SelectorAbsEtaRange(double absetamin, double absetamax);
966 Selector SelectorPhiRange(double phimin, double phimax);
967 Selector SelectorRapPhiRange(double rapmin, double rapmax, double phimin, double phimax);
968 Selector SelectorNHardest(unsigned int n);
969 Selector SelectorCircle(const double radius);
970 Selector SelectorDoughnut(const double radius_in, const double radius_out);
971 Selector SelectorStrip(const double half_width);
972 Selector SelectorRectangle(const double half_rap_width, const double half_phi_width);
973 Selector SelectorPtFractionMin(double fraction);
975 FJCORE_END_NAMESPACE // defined in fastjet/internal/base.hh
976 #endif // __FJCORE_SELECTOR_HH__
977 #ifndef __FJCORE_JETDEFINITION_HH__
978 #define __FJCORE_JETDEFINITION_HH__
979 #include<cassert>
980 #include<string>
981 #include<memory>
982 FJCORE_BEGIN_NAMESPACE // defined in fastjet/internal/base.hh
984 enum Strategy {
986  N2MHTLazy9 = -7,
987  N2MHTLazy25 = -6,
990  N2Tiled = -3,
991  N2PoorTiled = -2,
992  N2Plain = -1,
993  N3Dumb = 0,
994  Best = 1,
995  NlnN = 2,
996  NlnN3pi = 3,
997  NlnN4pi = 4,
1000  NlnNCam = 12, // 2piMultD
1001  BestFJ30 = 21,
1003 };
1015 };
1030 };
1031 class ClusterSequence;
1033 public:
1034  class Plugin;
1035  class Recombiner;
1036  JetDefinition(JetAlgorithm jet_algorithm_in,
1037  double R_in,
1038  RecombinationScheme recomb_scheme_in = E_scheme,
1039  Strategy strategy_in = Best) {
1040  *this = JetDefinition(jet_algorithm_in, R_in, recomb_scheme_in, strategy_in, 1);
1041  }
1042  JetDefinition(JetAlgorithm jet_algorithm_in,
1043  RecombinationScheme recomb_scheme_in = E_scheme,
1044  Strategy strategy_in = Best) {
1045  double dummyR = 0.0;
1046  *this = JetDefinition(jet_algorithm_in, dummyR, recomb_scheme_in, strategy_in, 0);
1047  }
1048  JetDefinition(JetAlgorithm jet_algorithm_in,
1049  double R_in,
1050  double xtra_param_in,
1051  RecombinationScheme recomb_scheme_in = E_scheme,
1052  Strategy strategy_in = Best) {
1053  *this = JetDefinition(jet_algorithm_in, R_in, recomb_scheme_in, strategy_in, 2);
1054  set_extra_param(xtra_param_in);
1055  }
1056  JetDefinition(JetAlgorithm jet_algorithm_in,
1057  double R_in,
1058  const Recombiner * recombiner_in,
1059  Strategy strategy_in = Best) {
1060  *this = JetDefinition(jet_algorithm_in, R_in, external_scheme, strategy_in);
1061  _recombiner = recombiner_in;
1062  }
1063  JetDefinition(JetAlgorithm jet_algorithm_in,
1064  const Recombiner * recombiner_in,
1065  Strategy strategy_in = Best) {
1066  *this = JetDefinition(jet_algorithm_in, external_scheme, strategy_in);
1067  _recombiner = recombiner_in;
1068  }
1069  JetDefinition(JetAlgorithm jet_algorithm_in,
1070  double R_in,
1071  double xtra_param_in,
1072  const Recombiner * recombiner_in,
1073  Strategy strategy_in = Best) {
1074  *this = JetDefinition(jet_algorithm_in, R_in, xtra_param_in, external_scheme, strategy_in);
1075  _recombiner = recombiner_in;
1076  }
1078  *this = JetDefinition(undefined_jet_algorithm, 1.0);
1079  }
1080  JetDefinition(const Plugin * plugin_in) {
1081  _plugin = plugin_in;
1083  _Rparam = _plugin->R();
1084  _extra_param = 0.0; // a dummy value to keep static code checkers happy
1087  }
1088  JetDefinition(JetAlgorithm jet_algorithm_in,
1089  double R_in,
1090  RecombinationScheme recomb_scheme_in,
1091  Strategy strategy_in,
1092  int nparameters_in);
1093  FJCORE_DEPRECATED_MSG("This argument ordering is deprecated. Use JetDefinition(alg, R, strategy, scheme[, n_parameters]) instead")
1094  JetDefinition(JetAlgorithm jet_algorithm_in,
1095  double R_in,
1096  Strategy strategy_in,
1097  RecombinationScheme recomb_scheme_in = E_scheme,
1098  int nparameters_in = 1){
1099  (*this) = JetDefinition(jet_algorithm_in,R_in,recomb_scheme_in,strategy_in,nparameters_in);
1100  }
1101  template <class L>
1102  std::vector<PseudoJet> operator()(const std::vector<L> & particles) const;
1103  static const double max_allowable_R; //= 1000.0;
1105  void set_recombiner(const Recombiner * recomb) {
1107  _recombiner = recomb;
1109  }
1110  void set_recombiner(const JetDefinition &other_jet_def);
1112  const Plugin * plugin() const {return _plugin;};
1116  double R () const {return _Rparam ;}
1117  double extra_param () const {return _extra_param ;}
1118  Strategy strategy () const {return _strategy ;}
1120  return _default_recombiner.scheme();}
1123  void set_extra_param(double xtra_param) {_extra_param = xtra_param;}
1124  const Recombiner * recombiner() const {
1125  return _recombiner == 0 ? & _default_recombiner : _recombiner;}
1126  bool has_same_recombiner(const JetDefinition &other_jd) const;
1127  bool is_spherical() const;
1128  std::string description() const;
1130  static std::string algorithm_description(const JetAlgorithm jet_alg);
1131  static unsigned int n_parameters_for_algorithm(const JetAlgorithm jet_alg);
1132 public:
1133  class Recombiner {
1134  public:
1135  virtual std::string description() const = 0;
1136  virtual void recombine(const PseudoJet & pa, const PseudoJet & pb,
1137  PseudoJet & pab) const = 0;
1138  virtual void preprocess(PseudoJet & ) const {};
1139  virtual ~Recombiner() {};
1140  inline void plus_equal(PseudoJet & pa, const PseudoJet & pb) const {
1141  PseudoJet pres;
1142  recombine(pa,pb,pres);
1143  pa = pres;
1144  }
1145  };
1147  public:
1149  _recomb_scheme(recomb_scheme) {}
1150  virtual std::string description() const FJCORE_OVERRIDE;
1151  virtual void recombine(const PseudoJet & pa, const PseudoJet & pb,
1152  PseudoJet & pab) const FJCORE_OVERRIDE;
1153  virtual void preprocess(PseudoJet & p) const FJCORE_OVERRIDE;
1155  private:
1157  };
1158  class Plugin{
1159  public:
1160  virtual std::string description() const = 0;
1161  virtual void run_clustering(ClusterSequence &) const = 0;
1162  virtual double R() const = 0;
1163  virtual bool supports_ghosted_passive_areas() const {return false;}
1164  virtual void set_ghost_separation_scale(double scale) const;
1165  virtual double ghost_separation_scale() const {return 0.0;}
1166  virtual bool exclusive_sequence_meaningful() const {return false;}
1167  virtual bool is_spherical() const {return false;}
1168  virtual ~Plugin() {};
1169  };
1170 private:
1172  double _Rparam;
1173  double _extra_param ;
1175  const Plugin * _plugin;
1180 };
1181 PseudoJet join(const std::vector<PseudoJet> & pieces, const JetDefinition::Recombiner & recombiner);
1182 PseudoJet join(const PseudoJet & j1,
1183  const JetDefinition::Recombiner & recombiner);
1184 PseudoJet join(const PseudoJet & j1, const PseudoJet & j2,
1185  const JetDefinition::Recombiner & recombiner);
1186 PseudoJet join(const PseudoJet & j1, const PseudoJet & j2, const PseudoJet & j3,
1187  const JetDefinition::Recombiner & recombiner);
1188 PseudoJet join(const PseudoJet & j1, const PseudoJet & j2, const PseudoJet & j3, const PseudoJet & j4,
1189  const JetDefinition::Recombiner & recombiner);
1191 #endif // __FJCORE_JETDEFINITION_HH__
1192 #ifndef __FJCORE_COMPOSITEJET_STRUCTURE_HH__
1193 #define __FJCORE_COMPOSITEJET_STRUCTURE_HH__
1194 FJCORE_BEGIN_NAMESPACE // defined in fastjet/internal/base.hh
1196 public:
1198  CompositeJetStructure(const std::vector<PseudoJet> & initial_pieces,
1199  const JetDefinition::Recombiner * recombiner = 0);
1202  };
1203  virtual std::string description() const FJCORE_OVERRIDE;
1204  virtual bool has_constituents() const FJCORE_OVERRIDE;
1205  virtual std::vector<PseudoJet> constituents(const PseudoJet &jet) const FJCORE_OVERRIDE;
1206  virtual bool has_pieces(const PseudoJet & /*jet*/) const FJCORE_OVERRIDE {return true;}
1207  virtual std::vector<PseudoJet> pieces(const PseudoJet &jet) const FJCORE_OVERRIDE;
1208 protected:
1209  std::vector<PseudoJet> _pieces;
1211 };
1212 template<typename T> PseudoJet join(const std::vector<PseudoJet> & pieces){
1213  PseudoJet result(0.0,0.0,0.0,0.0);
1214  for (unsigned int i=0; i<pieces.size(); i++){
1215  const PseudoJet it = pieces[i];
1216  result += it;
1217  }
1218  T *cj_struct = new T(pieces);
1220  return result;
1221 }
1222 template<typename T> PseudoJet join(const PseudoJet & j1){
1223  return join<T>(std::vector<PseudoJet>(1,j1));
1224 }
1225 template<typename T> PseudoJet join(const PseudoJet & j1, const PseudoJet & j2){
1226  std::vector<PseudoJet> pieces;
1227  pieces.push_back(j1);
1228  pieces.push_back(j2);
1229  return join<T>(pieces);
1230 }
1231 template<typename T> PseudoJet join(const PseudoJet & j1, const PseudoJet & j2,
1232  const PseudoJet & j3){
1233  std::vector<PseudoJet> pieces;
1234  pieces.push_back(j1);
1235  pieces.push_back(j2);
1236  pieces.push_back(j3);
1237  return join<T>(pieces);
1238 }
1239 template<typename T> PseudoJet join(const PseudoJet & j1, const PseudoJet & j2,
1240  const PseudoJet & j3, const PseudoJet & j4){
1241  std::vector<PseudoJet> pieces;
1242  pieces.push_back(j1);
1243  pieces.push_back(j2);
1244  pieces.push_back(j3);
1245  pieces.push_back(j4);
1246  return join<T>(pieces);
1247 }
1248 template<typename T> PseudoJet join(const std::vector<PseudoJet> & pieces,
1249  const JetDefinition::Recombiner & recombiner){
1250  PseudoJet result;
1251  if (pieces.size()>0){
1252  result = pieces[0];
1253  for (unsigned int i=1; i<pieces.size(); i++){
1254  recombiner.plus_equal(result, pieces[i]);
1255  }
1256  }
1257  T *cj_struct = new T(pieces, &recombiner);
1259  return result;
1260 }
1261 template<typename T> PseudoJet join(const PseudoJet & j1,
1262  const JetDefinition::Recombiner & recombiner){
1263  return join<T>(std::vector<PseudoJet>(1,j1), recombiner);
1264 }
1265 template<typename T> PseudoJet join(const PseudoJet & j1, const PseudoJet & j2,
1266  const JetDefinition::Recombiner & recombiner){
1267  std::vector<PseudoJet> pieces;
1268  pieces.reserve(2);
1269  pieces.push_back(j1);
1270  pieces.push_back(j2);
1271  return join<T>(pieces, recombiner);
1272 }
1273 template<typename T> PseudoJet join(const PseudoJet & j1, const PseudoJet & j2,
1274  const PseudoJet & j3,
1275  const JetDefinition::Recombiner & recombiner){
1276  std::vector<PseudoJet> pieces;
1277  pieces.reserve(3);
1278  pieces.push_back(j1);
1279  pieces.push_back(j2);
1280  pieces.push_back(j3);
1281  return join<T>(pieces, recombiner);
1282 }
1283 template<typename T> PseudoJet join(const PseudoJet & j1, const PseudoJet & j2,
1284  const PseudoJet & j3, const PseudoJet & j4,
1285  const JetDefinition::Recombiner & recombiner){
1286  std::vector<PseudoJet> pieces;
1287  pieces.reserve(4);
1288  pieces.push_back(j1);
1289  pieces.push_back(j2);
1290  pieces.push_back(j3);
1291  pieces.push_back(j4);
1292  return join<T>(pieces, recombiner);
1293 }
1294 FJCORE_END_NAMESPACE // defined in fastjet/internal/base.hh
1295 #endif // __FJCORE_MERGEDJET_STRUCTURE_HH__
1296 #ifndef __FJCORE_CLUSTER_SEQUENCE_STRUCTURE_HH__
1297 #define __FJCORE_CLUSTER_SEQUENCE_STRUCTURE_HH__
1298 #include <vector>
1299 FJCORE_BEGIN_NAMESPACE // defined in fastjet/internal/base.hh
1301 public:
1304  set_associated_cs(cs);
1305  };
1306  virtual ~ClusterSequenceStructure();
1308  return "PseudoJet with an associated ClusterSequence";
1309  }
1310  virtual bool has_associated_cluster_sequence() const FJCORE_OVERRIDE{ return true;}
1312  virtual bool has_valid_cluster_sequence() const FJCORE_OVERRIDE;
1313  virtual const ClusterSequence * validated_cs() const FJCORE_OVERRIDE;
1314  virtual void set_associated_cs(const ClusterSequence * new_cs){
1315  _associated_cs = new_cs;
1316  }
1317  virtual bool has_partner(const PseudoJet &reference, PseudoJet &partner) const FJCORE_OVERRIDE;
1318  virtual bool has_child(const PseudoJet &reference, PseudoJet &child) const FJCORE_OVERRIDE;
1319  virtual bool has_parents(const PseudoJet &reference, PseudoJet &parent1, PseudoJet &parent2) const FJCORE_OVERRIDE;
1320  virtual bool object_in_jet(const PseudoJet &reference, const PseudoJet &jet) const FJCORE_OVERRIDE;
1321  virtual bool has_constituents() const FJCORE_OVERRIDE;
1322  virtual std::vector<PseudoJet> constituents(const PseudoJet &reference) const FJCORE_OVERRIDE;
1323  virtual bool has_exclusive_subjets() const FJCORE_OVERRIDE;
1324  virtual std::vector<PseudoJet> exclusive_subjets(const PseudoJet &reference, const double & dcut) const FJCORE_OVERRIDE;
1325  virtual int n_exclusive_subjets(const PseudoJet &reference, const double & dcut) const FJCORE_OVERRIDE;
1326  virtual std::vector<PseudoJet> exclusive_subjets_up_to (const PseudoJet &reference, int nsub) const FJCORE_OVERRIDE;
1327  virtual double exclusive_subdmerge(const PseudoJet &reference, int nsub) const FJCORE_OVERRIDE;
1328  virtual double exclusive_subdmerge_max(const PseudoJet &reference, int nsub) const FJCORE_OVERRIDE;
1329  virtual bool has_pieces(const PseudoJet &reference) const FJCORE_OVERRIDE;
1330  virtual std::vector<PseudoJet> pieces(const PseudoJet &reference) const FJCORE_OVERRIDE;
1331 protected:
1333 };
1335 #endif // __FJCORE_CLUSTER_SEQUENCE_STRUCTURE_HH__
1336 #ifndef __FJCORE_CLUSTERSEQUENCE_HH__
1337 #define __FJCORE_CLUSTERSEQUENCE_HH__
1338 #include<vector>
1339 #include<map>
1340 #include<memory>
1341 #include<cassert>
1342 #include<iostream>
1343 #include<string>
1344 #include<set>
1345 #include<cmath> // needed to get double std::abs(double)
1346 FJCORE_BEGIN_NAMESPACE // defined in fastjet/internal/base.hh
1350  public:
1351  ClusterSequence () : _deletes_self_when_unused(false) {}
1352  template<class L> ClusterSequence (
1353  const std::vector<L> & pseudojets,
1354  const JetDefinition & jet_def,
1355  const bool & writeout_combinations = false);
1356  ClusterSequence (const ClusterSequence & cs) : _deletes_self_when_unused(false) {
1357  transfer_from_sequence(cs);
1358  }
1360  virtual ~ClusterSequence (); //{}
1361  std::vector<PseudoJet> inclusive_jets (const double ptmin = 0.0) const;
1362  int n_exclusive_jets (const double dcut) const;
1363  std::vector<PseudoJet> exclusive_jets (const double dcut) const;
1364  std::vector<PseudoJet> exclusive_jets (const int njets) const;
1365  std::vector<PseudoJet> exclusive_jets_up_to (const int njets) const;
1366  double exclusive_dmerge (const int njets) const;
1367  double exclusive_dmerge_max (const int njets) const;
1368  double exclusive_ymerge (int njets) const {return exclusive_dmerge(njets) / Q2();}
1369  double exclusive_ymerge_max (int njets) const {return exclusive_dmerge_max(njets)/Q2();}
1370  int n_exclusive_jets_ycut (double ycut) const {return n_exclusive_jets(ycut*Q2());}
1371  std::vector<PseudoJet> exclusive_jets_ycut (double ycut) const {
1372  int njets = n_exclusive_jets_ycut(ycut);
1373  return exclusive_jets(njets);
1374  }
1375  std::vector<PseudoJet> exclusive_subjets (const PseudoJet & jet,
1376  const double dcut) const;
1377  int n_exclusive_subjets(const PseudoJet & jet,
1378  const double dcut) const;
1379  std::vector<PseudoJet> exclusive_subjets (const PseudoJet & jet,
1380  int nsub) const;
1381  std::vector<PseudoJet> exclusive_subjets_up_to (const PseudoJet & jet,
1382  int nsub) const;
1383  double exclusive_subdmerge(const PseudoJet & jet, int nsub) const;
1384  double exclusive_subdmerge_max(const PseudoJet & jet, int nsub) const;
1385  double Q() const {return _Qtot;}
1386  double Q2() const {return _Qtot*_Qtot;}
1387  bool object_in_jet(const PseudoJet & object, const PseudoJet & jet) const;
1388  bool has_parents(const PseudoJet & jet, PseudoJet & parent1,
1389  PseudoJet & parent2) const;
1390  bool has_child(const PseudoJet & jet, PseudoJet & child) const;
1391  bool has_child(const PseudoJet & jet, const PseudoJet * & childp) const;
1392  bool has_partner(const PseudoJet & jet, PseudoJet & partner) const;
1393  std::vector<PseudoJet> constituents (const PseudoJet & jet) const;
1394  void print_jets_for_root(const std::vector<PseudoJet> & jets,
1395  std::ostream & ostr = std::cout) const;
1396  void print_jets_for_root(const std::vector<PseudoJet> & jets,
1397  const std::string & filename,
1398  const std::string & comment = "") const;
1399  void add_constituents (const PseudoJet & jet,
1400  std::vector<PseudoJet> & subjet_vector) const;
1401  inline Strategy strategy_used () const {return _strategy;}
1402  std::string strategy_string () const {return strategy_string(_strategy);}
1403  std::string strategy_string (Strategy strategy_in) const;
1404  const JetDefinition & jet_def() const {return _jet_def;}
1405  void delete_self_when_unused();
1406  bool will_delete_self_when_unused() const {return _deletes_self_when_unused;}
1407  void signal_imminent_self_deletion() const;
1408  double jet_scale_for_algorithm(const PseudoJet & jet) const;
1409  void plugin_record_ij_recombination(int jet_i, int jet_j, double dij,
1410  int & newjet_k) {
1411  assert(plugin_activated());
1412  _do_ij_recombination_step(jet_i, jet_j, dij, newjet_k);
1413  }
1414  void plugin_record_ij_recombination(int jet_i, int jet_j, double dij,
1415  const PseudoJet & newjet,
1416  int & newjet_k);
1417  void plugin_record_iB_recombination(int jet_i, double diB) {
1418  assert(plugin_activated());
1419  _do_iB_recombination_step(jet_i, diB);
1420  }
1421  class Extras {
1422  public:
1423  virtual ~Extras() {}
1424  virtual std::string description() const {return "This is a dummy extras class that contains no extra information! Derive from it if you want to use it to provide extra information from a plugin jet finder";}
1425  };
1426  inline void plugin_associate_extras(Extras * extras_in) {
1427  _extras.reset(extras_in);
1428  }
1429 #ifdef FJCORE_HAVE_AUTO_PTR_INTERFACE
1430  FJCORE_DEPRECATED_MSG("Please use ClusterSequence::plugin_associate_extras(Extras * extras_in)) instead")
1431  inline void plugin_associate_extras(std::auto_ptr<Extras> extras_in){
1432  _extras.reset(extras_in.release());
1433  }
1434 #endif
1435  inline bool plugin_activated() const {return _plugin_activated;}
1436  const Extras * extras() const {return _extras.get();}
1437  template<class GBJ> void plugin_simple_N2_cluster () {
1438  assert(plugin_activated());
1439  _simple_N2_cluster<GBJ>();
1440  }
1441 public:
1443  int parent1;
1444  int parent2;
1445  int child;
1446 
1447 
1448 
1450  double dij;
1451 
1453 
1454  };
1455  enum JetType {Invalid=-3, InexistentParent = -2, BeamJet = -1};
1456  const std::vector<PseudoJet> & jets() const;
1457  const std::vector<history_element> & history() const;
1458  unsigned int n_particles() const;
1459  std::vector<int> particle_jet_indices(const std::vector<PseudoJet> &) const;
1460  std::vector<int> unique_history_order() const;
1461  std::vector<PseudoJet> unclustered_particles() const;
1462  std::vector<PseudoJet> childless_pseudojets() const;
1463  bool contains(const PseudoJet & object) const;
1464  void transfer_from_sequence(const ClusterSequence & from_seq,
1465  const FunctionOfPseudoJet<PseudoJet> * action_on_jets = 0);
1467  return _structure_shared_ptr;
1468  }
1470  static void print_banner();
1471  static void set_fastjet_banner_stream(std::ostream * ostr) {_fastjet_banner_ostr = ostr;}
1472  static std::ostream * fastjet_banner_stream() {return _fastjet_banner_ostr;}
1473 private:
1474  static std::ostream * _fastjet_banner_ostr;
1475 protected:
1477  template<class L> void _transfer_input_jets(
1478  const std::vector<L> & pseudojets);
1479  void _initialise_and_run (const JetDefinition & jet_def,
1480  const bool & writeout_combinations);
1481  void _initialise_and_run_no_decant();
1482  void _decant_options(const JetDefinition & jet_def,
1483  const bool & writeout_combinations);
1484  void _decant_options_partial();
1485  void _fill_initial_history();
1486  void _do_ij_recombination_step(const int jet_i, const int jet_j,
1487  const double dij, int & newjet_k);
1488  void _do_iB_recombination_step(const int jet_i, const double diB);
1489  void _set_structure_shared_ptr(PseudoJet & j);
1490  void _update_structure_use_count();
1491  Strategy _best_strategy() const;
1492  class _Parabola {
1493  public:
1494  _Parabola(double a, double b, double c) : _a(a), _b(b), _c(c) {}
1495  inline double operator()(const double R) const {return _c*(_a*R*R + _b*R + 1);}
1496  private:
1497  double _a, _b, _c;
1498  };
1499  class _Line {
1500  public:
1501  _Line(double a, double b) : _a(a), _b(b) {}
1502  inline double operator()(const double R) const {return _a*R + _b;}
1503  private:
1504  double _a, _b;
1505  };
1506  std::vector<PseudoJet> _jets;
1507  std::vector<history_element> _history;
1508  void get_subhist_set(std::set<const history_element*> & subhist,
1509  const PseudoJet & jet, double dcut, int maxjet) const;
1512  double _Rparam, _R2, _invR2;
1513  double _Qtot;
1516  SharedPtr<PseudoJetStructureBase> _structure_shared_ptr; //< will actually be of type ClusterSequenceStructure
1517  int _structure_use_count_after_construction; //< info of use when CS handles its own memory
1519  private:
1521  SharedPtr<Extras> _extras; // things the plugin might want to add
1522  void _really_dumb_cluster ();
1523  void _delaunay_cluster ();
1524  template<class BJ> void _simple_N2_cluster ();
1525  void _tiled_N2_cluster ();
1526  void _faster_tiled_N2_cluster ();
1527  void _minheap_faster_tiled_N2_cluster();
1528  void _CP2DChan_cluster();
1529  void _CP2DChan_cluster_2pi2R ();
1530  void _CP2DChan_cluster_2piMultD ();
1531  void _CP2DChan_limited_cluster(double D);
1532  void _do_Cambridge_inclusive_jets();
1533  void _fast_NsqrtN_cluster();
1534  void _add_step_to_history( //const int step_number,
1535  const int parent1,
1536  const int parent2, const int jetp_index,
1537  const double dij);
1538  void _extract_tree_children(int pos, std::valarray<bool> &,
1539  const std::valarray<int> &, std::vector<int> &) const;
1540  void _extract_tree_parents (int pos, std::valarray<bool> &,
1541  const std::valarray<int> &, std::vector<int> &) const;
1542  typedef std::pair<int,int> TwoVertices;
1543  typedef std::pair<double,TwoVertices> DijEntry;
1544  typedef std::multimap<double,TwoVertices> DistMap;
1545  void _add_ktdistance_to_map(const int ii,
1546  DistMap & DijMap,
1547  const DynamicNearestNeighbours * DNN);
1548  static bool _first_time;
1551  struct BriefJet {
1552  double eta, phi, kt2, NN_dist;
1555  };
1556  class TiledJet {
1557  public:
1558  double eta, phi, kt2, NN_dist;
1560  int _jets_index, tile_index, diJ_posn;
1561  inline void label_minheap_update_needed() {diJ_posn = 1;}
1562  inline void label_minheap_update_done() {diJ_posn = 0;}
1563  inline bool minheap_update_needed() const {return diJ_posn==1;}
1564  };
1565  template <class J> void _bj_set_jetinfo( J * const jet,
1566  const int _jets_index) const;
1567  void _bj_remove_from_tiles( TiledJet * const jet) const;
1568  template <class J> double _bj_dist(const J * const jeta,
1569  const J * const jetb) const;
1570  template <class J> double _bj_diJ(const J * const jeta) const;
1571  template <class J> inline J * _bj_of_hindex(
1572  const int hist_index,
1573  J * const head, J * const tail)
1574  const {
1575  J * res;
1576  for(res = head; res<tail; res++) {
1577  if (_jets[res->_jets_index].cluster_hist_index() == hist_index) {break;}
1578  }
1579  return res;
1580  }
1581  template <class J> void _bj_set_NN_nocross(J * const jeta,
1582  J * const head, const J * const tail) const;
1583  template <class J> void _bj_set_NN_crosscheck(J * const jeta,
1584  J * const head, const J * const tail) const;
1585  static const int n_tile_neighbours = 9;
1586  struct Tile {
1587  Tile * begin_tiles[n_tile_neighbours];
1592  bool tagged;
1593  };
1594  std::vector<Tile> _tiles;
1595  double _tiles_eta_min, _tiles_eta_max;
1596  double _tile_size_eta, _tile_size_phi;
1597  int _n_tiles_phi,_tiles_ieta_min,_tiles_ieta_max;
1598  inline int _tile_index (int ieta, int iphi) const {
1599  return (ieta-_tiles_ieta_min)*_n_tiles_phi
1600  + (iphi+_n_tiles_phi) % _n_tiles_phi;
1601  }
1602  int _tile_index(const double eta, const double phi) const;
1603  void _tj_set_jetinfo ( TiledJet * const jet, const int _jets_index);
1604  void _bj_remove_from_tiles(TiledJet * const jet);
1605  void _initialise_tiles();
1606  void _print_tiles(TiledJet * briefjets ) const;
1607  void _add_neighbours_to_tile_union(const int tile_index,
1608  std::vector<int> & tile_union, int & n_near_tiles) const;
1609  void _add_untagged_neighbours_to_tile_union(const int tile_index,
1610  std::vector<int> & tile_union, int & n_near_tiles);
1611  struct EEBriefJet {
1612  double NN_dist; // obligatorily present
1613  double kt2; // obligatorily present == E^2 in general
1614  EEBriefJet * NN; // must be present too
1615  int _jets_index; // must also be present!
1616  double nx, ny, nz; // our internal storage for fast distance calcs
1617  };
1618  void _simple_N2_cluster_BriefJet();
1619  void _simple_N2_cluster_EEBriefJet();
1620 };
1622  const std::vector<L> & pseudojets) {
1623  _jets.reserve(pseudojets.size()*2);
1624  for (unsigned int i = 0; i < pseudojets.size(); i++) {
1625  _jets.push_back(pseudojets[i]);}
1626 }
1628  const std::vector<L> & pseudojets,
1629  const JetDefinition & jet_def_in,
1630  const bool & writeout_combinations) :
1631  _jet_def(jet_def_in), _writeout_combinations(writeout_combinations),
1632  _structure_shared_ptr(new ClusterSequenceStructure(this))
1633 {
1634  _transfer_input_jets(pseudojets);
1637 }
1638 inline const std::vector<PseudoJet> & ClusterSequence::jets () const {
1639  return _jets;
1640 }
1641 inline const std::vector<ClusterSequence::history_element> & ClusterSequence::history () const {
1642  return _history;
1643 }
1644 inline unsigned int ClusterSequence::n_particles() const {return _initial_n;}
1645 #ifndef __CINT__
1646 template<class L>
1647 std::vector<PseudoJet> JetDefinition::operator()(const std::vector<L> & particles) const {
1648  ClusterSequence * cs = new ClusterSequence(particles, *this);
1649  std::vector<PseudoJet> jets;
1650  if (is_spherical()) {
1651  jets = sorted_by_E(cs->inclusive_jets());
1652  } else {
1653  jets = sorted_by_pt(cs->inclusive_jets());
1654  }
1655  if (jets.size() != 0) {
1657  } else {
1658  delete cs;
1659  }
1660  return jets;
1661 }
1662 #endif // __CINT__
1663 template <class J> inline void ClusterSequence::_bj_set_jetinfo(
1664  J * const jetA, const int _jets_index) const {
1665  jetA->eta = _jets[_jets_index].rap();
1666  jetA->phi = _jets[_jets_index].phi_02pi();
1667  jetA->kt2 = jet_scale_for_algorithm(_jets[_jets_index]);
1668  jetA->_jets_index = _jets_index;
1669  jetA->NN_dist = _R2;
1670  jetA->NN = NULL;
1671 }
1672 template <class J> inline double ClusterSequence::_bj_dist(
1673  const J * const jetA, const J * const jetB) const {
1674 #ifndef FJCORE_NEW_DELTA_PHI
1675  double dphi = std::abs(jetA->phi - jetB->phi);
1676  double deta = (jetA->eta - jetB->eta);
1677  if (dphi > pi) {dphi = twopi - dphi;}
1678 #else
1679  double dphi = pi-std::abs(pi-std::abs(jetA->phi - jetB->phi));
1680  double deta = (jetA->eta - jetB->eta);
1681 #endif
1682  return dphi*dphi + deta*deta;
1683 }
1684 template <class J> inline double ClusterSequence::_bj_diJ(const J * const jet) const {
1685  double kt2 = jet->kt2;
1686  if (jet->NN != NULL) {if (jet->NN->kt2 < kt2) {kt2 = jet->NN->kt2;}}
1687  return jet->NN_dist * kt2;
1688 }
1689 template <class J> inline void ClusterSequence::_bj_set_NN_nocross(
1690  J * const jet, J * const head, const J * const tail) const {
1691  double NN_dist = _R2;
1692  J * NN = NULL;
1693  if (head < jet) {
1694  for (J * jetB = head; jetB != jet; jetB++) {
1695  double dist = _bj_dist(jet,jetB);
1696  if (dist < NN_dist) {
1697  NN_dist = dist;
1698  NN = jetB;
1699  }
1700  }
1701  }
1702  if (tail > jet) {
1703  for (J * jetB = jet+1; jetB != tail; jetB++) {
1704  double dist = _bj_dist(jet,jetB);
1705  if (dist < NN_dist) {
1706  NN_dist = dist;
1707  NN = jetB;
1708  }
1709  }
1710  }
1711  jet->NN = NN;
1712  jet->NN_dist = NN_dist;
1713 }
1714 template <class J> inline void ClusterSequence::_bj_set_NN_crosscheck(J * const jet,
1715  J * const head, const J * const tail) const {
1716  double NN_dist = _R2;
1717  J * NN = NULL;
1718  for (J * jetB = head; jetB != tail; jetB++) {
1719  double dist = _bj_dist(jet,jetB);
1720  if (dist < NN_dist) {
1721  NN_dist = dist;
1722  NN = jetB;
1723  }
1724  if (dist < jetB->NN_dist) {
1725  jetB->NN_dist = dist;
1726  jetB->NN = jet;
1727  }
1728  }
1729  jet->NN = NN;
1730  jet->NN_dist = NN_dist;
1731 }
1733 #endif // __FJCORE_CLUSTERSEQUENCE_HH__
1734 #ifndef __FJCORE_NNBASE_HH__
1735 #define __FJCORE_NNBASE_HH__
1736 FJCORE_BEGIN_NAMESPACE // defined in fastjet/internal/base.hh
1737 class _NoInfo {};
1738 template<class I> class NNInfo {
1739 public:
1740  NNInfo() : _info(NULL) {}
1741  NNInfo(I * info) : _info(info) {}
1742  template<class BJ> void init_jet(BJ * briefjet, const fjcore::PseudoJet & jet, int index) { briefjet->init(jet, index, _info);}
1743 private:
1745 };
1746 template<> class NNInfo<_NoInfo> {
1747 public:
1748  NNInfo() {}
1750  template<class BJ> void init_jet(BJ * briefjet, const fjcore::PseudoJet & jet, int index) { briefjet->init(jet, index);}
1751 };
1752 template<class I = _NoInfo> class NNBase : public NNInfo<I> {
1753 public:
1754  NNBase() {}
1755  NNBase(I * info) : NNInfo<I>(info) {}
1756  virtual void start(const std::vector<PseudoJet> & jets) = 0;
1757  virtual double dij_min(int & iA, int & iB) = 0;
1758  virtual void remove_jet(int iA) = 0;
1759  virtual void merge_jets(int iA, int iB, const PseudoJet & jet, int jet_index) = 0;
1760  virtual ~NNBase() {};
1761 };
1762 FJCORE_END_NAMESPACE // defined in fastjet/internal/base.hh
1763 #endif // __FJCORE_NNBASE_HH__
1764 #ifndef __FJCORE_NNH_HH__
1765 #define __FJCORE_NNH_HH__
1766 FJCORE_BEGIN_NAMESPACE // defined in fastjet/internal/base.hh
1767 template<class BJ, class I = _NoInfo> class NNH : public NNBase<I> {
1768 public:
1769  NNH(const std::vector<PseudoJet> & jets) : NNBase<I>() {start(jets);}
1770  NNH(const std::vector<PseudoJet> & jets, I * info) : NNBase<I>(info) {start(jets);}
1771  void start(const std::vector<PseudoJet> & jets);
1772  double dij_min(int & iA, int & iB);
1773  void remove_jet(int iA);
1774  void merge_jets(int iA, int iB, const PseudoJet & jet, int jet_index);
1775  ~NNH() {
1776  delete[] briefjets;
1777  }
1778 private:
1779  class NNBJ; // forward declaration
1780  void set_NN_crosscheck(NNBJ * jet, NNBJ * begin, NNBJ * end);
1781  void set_NN_nocross (NNBJ * jet, NNBJ * begin, NNBJ * end);
1784  int n;
1785  std::vector<NNBJ *> where_is;
1786  class NNBJ : public BJ {
1787  public:
1788  void init(const PseudoJet & jet, int index_in) {
1789  BJ::init(jet);
1790  other_init(index_in);
1791  }
1792  void init(const PseudoJet & jet, int index_in, I * info) {
1793  BJ::init(jet, info);
1794  other_init(index_in);
1795  }
1796  void other_init(int index_in) {
1797  _index = index_in;
1798  NN_dist = BJ::beam_distance();
1799  NN = NULL;
1800  }
1801  int index() const {return _index;}
1802  double NN_dist;
1804  private:
1805  int _index;
1806  };
1807 };
1808 template<class BJ, class I> void NNH<BJ,I>::start(const std::vector<PseudoJet> & jets) {
1809  n = jets.size();
1810  briefjets = new NNBJ[n];
1811  where_is.resize(2*n);
1812  NNBJ * jetA = briefjets;
1813  for (int i = 0; i< n; i++) {
1814  this->init_jet(jetA, jets[i], i);
1815  where_is[i] = jetA;
1816  jetA++; // move on to next entry of briefjets
1817  }
1818  tail = jetA; // a semaphore for the end of briefjets
1819  head = briefjets; // a nicer way of naming start
1820  for (jetA = head + 1; jetA != tail; jetA++) {
1821  set_NN_crosscheck(jetA, head, jetA);
1822  }
1823 }
1824 template<class BJ, class I> double NNH<BJ,I>::dij_min(int & iA, int & iB) {
1825  double diJ_min = briefjets[0].NN_dist;
1826  int diJ_min_jet = 0;
1827  for (int i = 1; i < n; i++) {
1828  if (briefjets[i].NN_dist < diJ_min) {
1829  diJ_min_jet = i;
1830  diJ_min = briefjets[i].NN_dist;
1831  }
1832  }
1833  NNBJ * jetA = & briefjets[diJ_min_jet];
1834  iA = jetA->index();
1835  iB = jetA->NN ? jetA->NN->index() : -1;
1836  return diJ_min;
1837 }
1838 template<class BJ, class I> void NNH<BJ,I>::remove_jet(int iA) {
1839  NNBJ * jetA = where_is[iA];
1840  tail--; n--;
1841  *jetA = *tail;
1842  where_is[jetA->index()] = jetA;
1843  for (NNBJ * jetI = head; jetI != tail; jetI++) {
1844  if (jetI->NN == jetA) set_NN_nocross(jetI, head, tail);
1845  if (jetI->NN == tail) {jetI->NN = jetA;}
1846  }
1847 }
1848 template<class BJ, class I> void NNH<BJ,I>::merge_jets(int iA, int iB,
1849  const PseudoJet & jet, int index) {
1850  NNBJ * jetA = where_is[iA];
1851  NNBJ * jetB = where_is[iB];
1852  if (jetA < jetB) std::swap(jetA,jetB);
1853  this->init_jet(jetB, jet, index);
1854  if (index >= int(where_is.size())) where_is.resize(2*index);
1855  where_is[jetB->index()] = jetB;
1856  tail--; n--;
1857  *jetA = *tail;
1858  where_is[jetA->index()] = jetA;
1859  for (NNBJ * jetI = head; jetI != tail; jetI++) {
1860  if (jetI->NN == jetA || jetI->NN == jetB) {
1861  set_NN_nocross(jetI, head, tail);
1862  }
1863  double dist = jetI->distance(jetB);
1864  if (dist < jetI->NN_dist) {
1865  if (jetI != jetB) {
1866  jetI->NN_dist = dist;
1867  jetI->NN = jetB;
1868  }
1869  }
1870  if (dist < jetB->NN_dist) {
1871  if (jetI != jetB) {
1872  jetB->NN_dist = dist;
1873  jetB->NN = jetI;
1874  }
1875  }
1876  if (jetI->NN == tail) {jetI->NN = jetA;}
1877  }
1878 }
1879 template <class BJ, class I> void NNH<BJ,I>::set_NN_crosscheck(NNBJ * jet,
1880  NNBJ * begin, NNBJ * end) {
1881  double NN_dist = jet->beam_distance();
1882  NNBJ * NN = NULL;
1883  for (NNBJ * jetB = begin; jetB != end; jetB++) {
1884  double dist = jet->distance(jetB);
1885  if (dist < NN_dist) {
1886  NN_dist = dist;
1887  NN = jetB;
1888  }
1889  if (dist < jetB->NN_dist) {
1890  jetB->NN_dist = dist;
1891  jetB->NN = jet;
1892  }
1893  }
1894  jet->NN = NN;
1895  jet->NN_dist = NN_dist;
1896 }
1897 template <class BJ, class I> void NNH<BJ,I>::set_NN_nocross(
1898  NNBJ * jet, NNBJ * begin, NNBJ * end) {
1899  double NN_dist = jet->beam_distance();
1900  NNBJ * NN = NULL;
1901  if (begin < jet) {
1902  for (NNBJ * jetB = begin; jetB != jet; jetB++) {
1903  double dist = jet->distance(jetB);
1904  if (dist < NN_dist) {
1905  NN_dist = dist;
1906  NN = jetB;
1907  }
1908  }
1909  }
1910  if (end > jet) {
1911  for (NNBJ * jetB = jet+1; jetB != end; jetB++) {
1912  double dist = jet->distance (jetB);
1913  if (dist < NN_dist) {
1914  NN_dist = dist;
1915  NN = jetB;
1916  }
1917  }
1918  }
1919  jet->NN = NN;
1920  jet->NN_dist = NN_dist;
1921 }
1922 FJCORE_END_NAMESPACE // defined in fastjet/internal/base.hh
1923 #endif // __FJCORE_NNH_HH__
1924 #endif