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