Line data Source code
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; ///< the value (only member carrying info)
102 : typedef T value_type; ///< a typedef for the type T
103 : typedef integral_type<T,_t> type; ///< a typedef for the whole structure
104 : };
105 : template<typename T, T _t>
106 : const T integral_type<T, _t>::value;
107 : typedef integral_type<bool, true> true_type; ///< the bool 'true' value promoted to a type
108 : typedef integral_type<bool, false> false_type; ///< the bool 'false' value promoted to a 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 0 : SharedPtr() : _ptr(NULL){}
242 0 : template<class Y> explicit SharedPtr(Y* ptr){
243 0 : _ptr = new __SharedCountingPtr(ptr);
244 0 : }
245 0 : SharedPtr(SharedPtr const & share) : _ptr(share._get_container()){
246 0 : if (_ptr!=NULL) ++(*_ptr);
247 0 : }
248 0 : ~SharedPtr(){
249 0 : if (_ptr==NULL) return;
250 0 : _decrease_count();
251 0 : }
252 : void reset(){
253 0 : SharedPtr().swap(*this);
254 0 : }
255 : template<class Y> void reset(Y * ptr){
256 0 : SharedPtr(ptr).swap(*this);
257 0 : }
258 : template<class Y> void reset(SharedPtr<Y> const & share){
259 0 : if (_ptr!=NULL){
260 0 : if (_ptr == share._get_container()) return;
261 0 : _decrease_count();
262 0 : }
263 0 : _ptr = share._get_container(); // Note: automatically set it to NULL if share is empty
264 0 : if (_ptr!=NULL) ++(*_ptr);
265 0 : }
266 : SharedPtr& operator=(SharedPtr const & share){
267 0 : reset(share);
268 0 : 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 0 : if (_ptr==NULL) return NULL;
276 0 : return _ptr->get(); // automatically returns NULL when out-of-scope
277 0 : }
278 : inline T& operator*() const{
279 : return *(_ptr->get());
280 : }
281 : inline T* operator->() const{
282 0 : if (_ptr==NULL) return NULL;
283 0 : return _ptr->get();
284 0 : }
285 : inline T* get() const{
286 0 : if (_ptr==NULL) return NULL;
287 0 : return _ptr->get();
288 0 : }
289 : inline bool unique() const{
290 0 : return (use_count()==1);
291 : }
292 : inline long use_count() const{
293 0 : if (_ptr==NULL) return 0;
294 0 : return _ptr->use_count(); // automatically returns NULL when out-of-scope
295 0 : }
296 : inline operator bool() const{
297 : return (get()!=NULL);
298 : }
299 : inline void swap(SharedPtr & share){
300 0 : __SharedCountingPtr* share_container = share._ptr;
301 0 : share._ptr = _ptr;
302 0 : _ptr = share_container;
303 0 : }
304 : void set_count(const long & count){
305 0 : if (_ptr==NULL) return;
306 0 : _ptr->set_count(count);
307 0 : }
308 : class __SharedCountingPtr{
309 : public:
310 : __SharedCountingPtr() : _ptr(NULL), _count(0){}
311 0 : template<class Y> explicit __SharedCountingPtr(Y* ptr) : _ptr(ptr), _count(1){}
312 0 : ~__SharedCountingPtr(){
313 0 : if (_ptr!=NULL){ delete _ptr;}
314 0 : }
315 0 : inline T* get() const {return _ptr;}
316 0 : inline long use_count() const {return _count;}
317 0 : inline long operator++(){return ++_count;}
318 : inline long operator--(){return --_count;}
319 : inline long operator++(int){return _count++;}
320 0 : inline long operator--(int){return _count--;}
321 : void set_count(const long & count){
322 0 : _count = count;
323 0 : }
324 : private:
325 : T *_ptr; ///< the pointer we're counting the references to
326 : long _count; ///< the number of references
327 : };
328 : private:
329 : inline __SharedCountingPtr* _get_container() const{
330 0 : return _ptr;
331 : }
332 : void _decrease_count(){
333 0 : (*_ptr)--;
334 0 : if (_ptr->use_count()==0)
335 0 : delete _ptr; // that automatically deletes the object itself
336 0 : }
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 0 : 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; ///< error message
380 : static bool _print_errors; ///< do we print anything?
381 : static bool _print_backtrace; ///< do we print the backtrace?
382 : static std::ostream * _default_ostr; ///< the output stream (cerr if not set)
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 12 : 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 0 : PseudoJetStructureBase(){};
425 0 : virtual ~PseudoJetStructureBase(){};
426 0 : virtual std::string description() const{ return "PseudoJet with an unknown structure"; }
427 0 : virtual bool has_associated_cluster_sequence() const { return false;}
428 : virtual const ClusterSequence* associated_cluster_sequence() const;
429 0 : 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 0 : virtual bool has_constituents() const {return false;}
436 : virtual std::vector<PseudoJet> constituents(const PseudoJet &reference) const;
437 0 : 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 0 : 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 0 : class PseudoJet {
462 : public:
463 0 : 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 0 : PseudoJet(bool /* dummy */) {}
467 0 : virtual ~PseudoJet(){};
468 0 : inline double E() const {return _E;}
469 0 : inline double e() const {return _E;} // like CLHEP
470 0 : inline double px() const {return _px;}
471 0 : inline double py() const {return _py;}
472 0 : inline double pz() const {return _pz;}
473 0 : inline double phi() const {return phi_02pi();}
474 : inline double phi_std() const {
475 0 : _ensure_valid_rap_phi();
476 0 : return _phi > pi ? _phi-twopi : _phi;}
477 : inline double phi_02pi() const {
478 0 : _ensure_valid_rap_phi();
479 0 : return _phi;
480 : }
481 : inline double rap() const {
482 0 : _ensure_valid_rap_phi();
483 0 : return _rap;
484 : }
485 : inline double rapidity() const {return rap();} // like CLHEP
486 : double pseudorapidity() const;
487 0 : double eta() const {return pseudorapidity();}
488 0 : inline double pt2() const {return _kt2;}
489 : inline double pt() const {return sqrt(_kt2);}
490 0 : inline double perp2() const {return _kt2;} // like CLHEP
491 0 : inline double perp() const {return sqrt(_kt2);} // like CLHEP
492 0 : inline double kt2() const {return _kt2;} // for bkwds compatibility
493 0 : 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 0 : 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 0 : 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 0 : 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 0 : reset_momentum_PtYPhiM(pt_in, y_in, phi_in, m_in);
537 0 : _reset_indices();
538 0 : }
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 0 : inline int user_index() const {return _user_index;}
548 0 : inline void set_user_index(const int index) {_user_index = index;}
549 : class UserInfoBase{
550 : public:
551 : UserInfoBase(){};
552 : virtual ~UserInfoBase(){};
553 : };
554 0 : 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 0 : if (!_user_info()) return NULL;
575 0 : return _user_info.get();
576 0 : }
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 0 : inline int cluster_hist_index() const {return _cluster_hist_index;}
623 0 : 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 0 : if (_phi == pseudojet_invalid_phi) _set_rap_phi();
640 0 : }
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 0 : inline IndexedSortHelper (const std::vector<double> * reference_values) {
668 0 : _ref_values = reference_values;
669 0 : };
670 : inline int operator() (const int & i1, const int & i2) const {
671 0 : 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 0 : set_cluster_hist_index(-1);
681 0 : set_user_index(-1);
682 0 : _structure.reset();
683 0 : _user_info.reset();
684 0 : }
685 : inline double PseudoJet::m() const {
686 0 : double mm = m2();
687 0 : 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 0 : _px = px_in;
691 0 : _py = py_in;
692 0 : _pz = pz_in;
693 0 : _E = E_in;
694 0 : _finish_init();
695 0 : _reset_indices();
696 0 : }
697 : inline void PseudoJet::reset_momentum(double px_in, double py_in, double pz_in, double E_in) {
698 0 : _px = px_in;
699 0 : _py = py_in;
700 0 : _pz = pz_in;
701 0 : _E = E_in;
702 0 : _finish_init();
703 0 : }
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 0 : std::vector<TOut> res(pjs.size());
749 0 : for (unsigned int i=0; i<pjs.size(); i++)
750 0 : res[i] = result(pjs[i]);
751 : return res;
752 0 : }
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 0 : class SelectorWorker {
763 : public:
764 0 : virtual ~SelectorWorker() {}
765 : virtual bool pass(const PseudoJet & jet) const = 0;
766 : virtual void terminator(std::vector<const PseudoJet *> & jets) const {
767 0 : for (unsigned i = 0; i < jets.size(); i++) {
768 0 : if (jets[i] && !pass(*jets[i])) jets[i] = NULL;
769 : }
770 0 : }
771 0 : virtual bool applies_jet_by_jet() const {return true;}
772 0 : virtual std::string description() const {return "missing description";}
773 0 : virtual bool takes_reference() const { return false;}
774 : virtual void set_reference(const PseudoJet & /*reference*/){
775 0 : throw Error("set_reference(...) cannot be used for a selector worker that does not take a reference");
776 0 : }
777 : virtual SelectorWorker* copy(){
778 0 : throw Error("this SelectorWorker has nothing to copy");
779 0 : }
780 : virtual void get_rapidity_extent(double & rapmin, double & rapmax) const {
781 0 : rapmax = std::numeric_limits<double>::infinity();
782 0 : rapmin = -rapmax;
783 0 : }
784 0 : virtual bool is_geometric() const { return false;}
785 : virtual bool has_finite_area() const;
786 0 : virtual bool has_known_area() const { return false;}
787 : virtual double known_area() const{
788 0 : throw Error("this selector has no computable area");
789 0 : }
790 : };
791 0 : class Selector{
792 : public:
793 : Selector() {}
794 0 : Selector(SelectorWorker * worker_in) {_worker.reset(worker_in);}
795 0 : virtual ~Selector(){}
796 : bool pass(const PseudoJet & jet) const {
797 0 : if (!validated_worker()->applies_jet_by_jet()) {
798 0 : throw Error("Cannot apply this selector to an individual jet");
799 : }
800 0 : return _worker->pass(jet);
801 0 : }
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 0 : 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 0 : validated_worker()->terminator(jets);
815 0 : }
816 : void get_rapidity_extent(double &rapmin, double &rapmax) const {
817 0 : return validated_worker()->get_rapidity_extent(rapmin, rapmax);
818 : }
819 : std::string description() const {
820 0 : return validated_worker()->description();
821 : }
822 : bool is_geometric() const{
823 0 : return validated_worker()->is_geometric();
824 : }
825 : bool has_finite_area() const{
826 : return validated_worker()->has_finite_area();
827 : }
828 0 : const SharedPtr<SelectorWorker> & worker() const {return _worker;}
829 : const SelectorWorker* validated_worker() const {
830 0 : const SelectorWorker* worker_ptr = _worker.get();
831 0 : if (worker_ptr == 0) throw InvalidWorker();
832 0 : return worker_ptr;
833 0 : }
834 : bool takes_reference() const {
835 0 : return validated_worker()->takes_reference();
836 : }
837 : const Selector & set_reference(const PseudoJet &reference){
838 0 : if (! validated_worker()->takes_reference()){
839 0 : return *this;
840 : }
841 0 : _copy_worker_if_needed();
842 0 : _worker->set_reference(reference);
843 0 : return *this;
844 0 : }
845 0 : class InvalidWorker : public Error {
846 : public:
847 0 : 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 0 : if (_worker.unique()) return;
858 0 : _worker.reset(_worker->copy());
859 0 : }
860 : private:
861 : SharedPtr<SelectorWorker> _worker; ///< the underlying 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); ///< select objects with pt >= ptmin
869 : Selector SelectorPtMax(double ptmax); ///< select objects with pt <= ptmax
870 : Selector SelectorPtRange(double ptmin, double ptmax); ///< select objects with ptmin <= pt <= ptmax
871 : Selector SelectorEtMin(double Etmin); ///< select objects with Et >= Etmin
872 : Selector SelectorEtMax(double Etmax); ///< select objects with Et <= Etmax
873 : Selector SelectorEtRange(double Etmin, double Etmax); ///< select objects with Etmin <= Et <= Etmax
874 : Selector SelectorEMin(double Emin); ///< select objects with E >= Emin
875 : Selector SelectorEMax(double Emax); ///< select objects with E <= Emax
876 : Selector SelectorERange(double Emin, double Emax); ///< select objects with Emin <= E <= Emax
877 : Selector SelectorMassMin(double Mmin); ///< select objects with Mass >= Mmin
878 : Selector SelectorMassMax(double Mmax); ///< select objects with Mass <= Mmax
879 : Selector SelectorMassRange(double Mmin, double Mmax); ///< select objects with Mmin <= Mass <= Mmax
880 : Selector SelectorRapMin(double rapmin); ///< select objects with rap >= rapmin
881 : Selector SelectorRapMax(double rapmax); ///< select objects with rap <= rapmax
882 : Selector SelectorRapRange(double rapmin, double rapmax); ///< select objects with rapmin <= rap <= rapmax
883 : Selector SelectorAbsRapMin(double absrapmin); ///< select objects with |rap| >= absrapmin
884 : Selector SelectorAbsRapMax(double absrapmax); ///< select objects with |rap| <= absrapmax
885 : Selector SelectorAbsRapRange(double absrapmin, double absrapmax); ///< select objects with absrapmin <= |rap| <= absrapmax
886 : Selector SelectorEtaMin(double etamin); ///< select objects with eta >= etamin
887 : Selector SelectorEtaMax(double etamax); ///< select objects with eta <= etamax
888 : Selector SelectorEtaRange(double etamin, double etamax); ///< select objects with etamin <= eta <= etamax
889 : Selector SelectorAbsEtaMin(double absetamin); ///< select objects with |eta| >= absetamin
890 : Selector SelectorAbsEtaMax(double absetamax); ///< select objects with |eta| <= absetamax
891 : Selector SelectorAbsEtaRange(double absetamin, double absetamax); ///< select objects with absetamin <= |eta| <= absetamax
892 : Selector SelectorPhiRange(double phimin, double phimax); ///< select objects with phimin <= phi <= 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 0 : class JetDefinition {
952 : public:
953 : class Plugin;
954 : class Recombiner;
955 0 : JetDefinition(JetAlgorithm jet_algorithm_in,
956 : double R_in,
957 : RecombinationScheme recomb_scheme_in = E_scheme,
958 0 : Strategy strategy_in = Best) {
959 0 : *this = JetDefinition(jet_algorithm_in, R_in, strategy_in, recomb_scheme_in, 1);
960 0 : }
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 0 : const Plugin * plugin() const {return _plugin;};
1021 : void delete_plugin_when_unused();
1022 0 : JetAlgorithm jet_algorithm () const {return _jet_algorithm ;}
1023 : JetAlgorithm jet_finder () const {return _jet_algorithm ;}
1024 0 : double R () const {return _Rparam ;}
1025 0 : double extra_param () const {return _extra_param ;}
1026 0 : Strategy strategy () const {return _strategy ;}
1027 : RecombinationScheme recombination_scheme() const {
1028 0 : 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 0 : void set_extra_param(double xtra_param) {_extra_param = xtra_param;}
1032 : const Recombiner * recombiner() const {
1033 0 : return _recombiner == 0 ? & _default_recombiner : _recombiner;}
1034 : bool has_same_recombiner(const JetDefinition &other_jd) const;
1035 : std::string description() const;
1036 : public:
1037 0 : 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 0 : virtual void preprocess(PseudoJet & ) const {};
1043 0 : virtual ~Recombiner() {};
1044 : inline void plus_equal(PseudoJet & pa, const PseudoJet & pb) const {
1045 0 : PseudoJet pres;
1046 0 : recombine(pa,pb,pres);
1047 0 : pa = pres;
1048 0 : }
1049 : };
1050 0 : class DefaultRecombiner : public Recombiner {
1051 : public:
1052 0 : DefaultRecombiner(RecombinationScheme recomb_scheme = E_scheme) :
1053 0 : _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 0 : 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 0 : virtual bool supports_ghosted_passive_areas() const {return false;}
1068 : virtual void set_ghost_separation_scale(double scale) const;
1069 0 : virtual double ghost_separation_scale() const {return 0.0;}
1070 0 : virtual bool exclusive_sequence_meaningful() const {return false;}
1071 0 : virtual ~Plugin() {};
1072 : };
1073 : private:
1074 : JetAlgorithm _jet_algorithm;
1075 : double _Rparam;
1076 : double _extra_param ; ///< parameter whose meaning varies according to context
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 0 : virtual ~CompositeJetStructure(){
1104 0 : if (_area_4vector_ptr) delete _area_4vector_ptr;
1105 0 : };
1106 : virtual std::string description() const;
1107 : virtual bool has_constituents() const;
1108 : virtual std::vector<PseudoJet> constituents(const PseudoJet &jet) const;
1109 0 : 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; ///< the pieces building the jet
1113 : PseudoJet * _area_4vector_ptr; ///< pointer to the 4-vector jet area
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 0 : ClusterSequenceStructure(const ClusterSequence *cs){
1204 0 : set_associated_cs(cs);
1205 0 : };
1206 : virtual ~ClusterSequenceStructure();
1207 0 : virtual std::string description() const{ return "PseudoJet with an associated ClusterSequence"; }
1208 0 : 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 0 : _associated_cs = new_cs;
1214 0 : }
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 0 : std::string strategy_string () const {return strategy_string(_strategy);}
1300 : std::string strategy_string (Strategy strategy_in) const;
1301 0 : const JetDefinition & jet_def() const {return _jet_def;}
1302 : void delete_self_when_unused();
1303 0 : 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 0 : assert(plugin_activated());
1309 0 : _do_ij_recombination_step(jet_i, jet_j, dij, newjet_k);
1310 0 : }
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 0 : 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; /// index in _history where first parent of this jet
1335 : int parent2; /// index in _history where second parent of this jet
1336 : int child; /// index in _history where the current jet is
1337 : /// recombined with another jet to form its child. It
1338 : /// is Invalid if this jet does not further
1339 : /// recombine.
1340 : int jetp_index; /// index in the _jets vector where we will find the
1341 : double dij; /// the distance corresponding to the recombination
1342 : /// at this stage of the clustering.
1343 : double max_dij_so_far; /// the largest recombination distance seen
1344 : /// so far in the clustering history.
1345 : };
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 0 : inline void label_minheap_update_needed() {diJ_posn = 1;}
1437 0 : inline void label_minheap_update_done() {diJ_posn = 0;}
1438 0 : 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 0 : return (ieta-_tiles_ieta_min)*_n_tiles_phi
1475 0 : + (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 0 : _jets.reserve(pseudojets.size()*2);
1499 0 : for (unsigned int i = 0; i < pseudojets.size(); i++) {
1500 0 : _jets.push_back(pseudojets[i]);}
1501 0 : }
1502 0 : template<class L> ClusterSequence::ClusterSequence (
1503 : const std::vector<L> & pseudojets,
1504 : const JetDefinition & jet_def_in,
1505 : const bool & writeout_combinations) :
1506 0 : _jet_def(jet_def_in), _writeout_combinations(writeout_combinations),
1507 0 : _structure_shared_ptr(new ClusterSequenceStructure(this))
1508 0 : {
1509 0 : _transfer_input_jets(pseudojets);
1510 0 : _decant_options_partial();
1511 0 : _initialise_and_run_no_decant();
1512 0 : }
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 0 : return _history;
1518 : }
1519 0 : 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 0 : jetA->eta = _jets[_jets_index].rap();
1523 0 : jetA->phi = _jets[_jets_index].phi_02pi();
1524 0 : jetA->kt2 = jet_scale_for_algorithm(_jets[_jets_index]);
1525 0 : jetA->_jets_index = _jets_index;
1526 0 : jetA->NN_dist = _R2;
1527 0 : jetA->NN = NULL;
1528 0 : }
1529 : template <class J> inline double ClusterSequence::_bj_dist(
1530 : const J * const jetA, const J * const jetB) const {
1531 0 : double dphi = std::abs(jetA->phi - jetB->phi);
1532 0 : double deta = (jetA->eta - jetB->eta);
1533 0 : if (dphi > pi) {dphi = twopi - dphi;}
1534 0 : return dphi*dphi + deta*deta;
1535 : }
1536 : template <class J> inline double ClusterSequence::_bj_diJ(const J * const jet) const {
1537 0 : double kt2 = jet->kt2;
1538 0 : if (jet->NN != NULL) {if (jet->NN->kt2 < kt2) {kt2 = jet->NN->kt2;}}
1539 0 : 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 0 : double NN_dist = _R2;
1544 : J * NN = NULL;
1545 0 : if (head < jet) {
1546 0 : for (J * jetB = head; jetB != jet; jetB++) {
1547 0 : double dist = _bj_dist(jet,jetB);
1548 0 : if (dist < NN_dist) {
1549 : NN_dist = dist;
1550 : NN = jetB;
1551 0 : }
1552 : }
1553 0 : }
1554 0 : if (tail > jet) {
1555 0 : for (J * jetB = jet+1; jetB != tail; jetB++) {
1556 0 : double dist = _bj_dist(jet,jetB);
1557 0 : if (dist < NN_dist) {
1558 : NN_dist = dist;
1559 : NN = jetB;
1560 0 : }
1561 : }
1562 0 : }
1563 0 : jet->NN = NN;
1564 0 : jet->NN_dist = NN_dist;
1565 0 : }
1566 : template <class J> inline void ClusterSequence::_bj_set_NN_crosscheck(J * const jet,
1567 : J * const head, const J * const tail) const {
1568 0 : double NN_dist = _R2;
1569 : J * NN = NULL;
1570 0 : for (J * jetB = head; jetB != tail; jetB++) {
1571 0 : double dist = _bj_dist(jet,jetB);
1572 0 : if (dist < NN_dist) {
1573 : NN_dist = dist;
1574 : NN = jetB;
1575 0 : }
1576 0 : if (dist < jetB->NN_dist) {
1577 0 : jetB->NN_dist = dist;
1578 0 : jetB->NN = jet;
1579 0 : }
1580 : }
1581 0 : jet->NN = NN;
1582 0 : jet->NN_dist = NN_dist;
1583 0 : }
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
|
