FastJet 3.4.0
ClusterSequence.hh
1#ifndef __FASTJET_CLUSTERSEQUENCE_HH__
2#define __FASTJET_CLUSTERSEQUENCE_HH__
3
4//FJSTARTHEADER
5// $Id$
6//
7// Copyright (c) 2005-2021, Matteo Cacciari, Gavin P. Salam and Gregory Soyez
8//
9//----------------------------------------------------------------------
10// This file is part of FastJet.
11//
12// FastJet is free software; you can redistribute it and/or modify
13// it under the terms of the GNU General Public License as published by
14// the Free Software Foundation; either version 2 of the License, or
15// (at your option) any later version.
16//
17// The algorithms that underlie FastJet have required considerable
18// development. They are described in the original FastJet paper,
19// hep-ph/0512210 and in the manual, arXiv:1111.6097. If you use
20// FastJet as part of work towards a scientific publication, please
21// quote the version you use and include a citation to the manual and
22// optionally also to hep-ph/0512210.
23//
24// FastJet is distributed in the hope that it will be useful,
25// but WITHOUT ANY WARRANTY; without even the implied warranty of
26// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
27// GNU General Public License for more details.
28//
29// You should have received a copy of the GNU General Public License
30// along with FastJet. If not, see <http://www.gnu.org/licenses/>.
31//----------------------------------------------------------------------
32//FJENDHEADER
33
34
35#include<vector>
36#include<map>
37#include "fastjet/PseudoJet.hh"
38#include<memory>
39#include<cassert>
40#include<iostream>
41#include<string>
42#include<set>
43#include<cmath> // needed to get double std::abs(double)
44#include "fastjet/Error.hh"
45#include "fastjet/JetDefinition.hh"
46#include "fastjet/SharedPtr.hh"
47#include "fastjet/LimitedWarning.hh"
48#include "fastjet/FunctionOfPseudoJet.hh"
49#include "fastjet/ClusterSequenceStructure.hh"
50#include "fastjet/internal/thread_safety_helpers.hh" // helpers to write code w&wo thread-safety
51#include "fastjet/internal/deprecated.hh"
52
53FASTJET_BEGIN_NAMESPACE // defined in fastjet/internal/base.hh
54
55
56// forward declarations
57class ClusterSequenceStructure;
58class DynamicNearestNeighbours;
59
60/// @ingroup basic_classes
61/// \class ClusterSequence
62/// deals with clustering
64
65
66 public:
67
68 /// default constructor
69 ClusterSequence () : _deletes_self_when_unused(false) {}
70
71 /// create a ClusterSequence, starting from the supplied set
72 /// of PseudoJets and clustering them with jet definition specified
73 /// by jet_def (which also specifies the clustering strategy)
74 template<class L> ClusterSequence (
75 const std::vector<L> & pseudojets,
76 const JetDefinition & jet_def,
77 const bool & writeout_combinations = false);
78
79 /// copy constructor for a ClusterSequence
80 ClusterSequence (const ClusterSequence & cs) : _deletes_self_when_unused(false) {
81 transfer_from_sequence(cs);
82 }
83
84 /// explicit assignment operator for a ClusterSequence
85 ClusterSequence & operator=(const ClusterSequence & cs);
86
87 // virtual ClusterSequence destructor, in case any derived class
88 // thinks of needing a destructor at some point
89 virtual ~ClusterSequence (); //{}
90
91 // NB: in the routines that follow, for extracting lists of jets, a
92 // list structure might be more efficient, if sometimes a little
93 // more awkward to use (at least for old fortran hands).
94
95 /// return a vector of all jets (in the sense of the inclusive
96 /// algorithm) with pt >= ptmin. Time taken should be of the order
97 /// of the number of jets returned.
98 std::vector<PseudoJet> inclusive_jets (const double ptmin = 0.0) const;
99
100 /// return the number of jets (in the sense of the exclusive
101 /// algorithm) that would be obtained when running the algorithm
102 /// with the given dcut.
103 int n_exclusive_jets (const double dcut) const;
104
105 /// return a vector of all jets (in the sense of the exclusive
106 /// algorithm) that would be obtained when running the algorithm
107 /// with the given dcut.
108 std::vector<PseudoJet> exclusive_jets (const double dcut) const;
109
110 /// return a vector of all jets when the event is clustered (in the
111 /// exclusive sense) to exactly njets.
112 ///
113 /// If there are fewer than njets particles in the ClusterSequence
114 /// an error is thrown
115 std::vector<PseudoJet> exclusive_jets (const int njets) const;
116
117 /// return a vector of all jets when the event is clustered (in the
118 /// exclusive sense) to exactly njets.
119 ///
120 /// If there are fewer than njets particles in the ClusterSequence
121 /// the function just returns however many particles there were.
122 std::vector<PseudoJet> exclusive_jets_up_to (const int njets) const;
123
124 /// return the dmin corresponding to the recombination that went
125 /// from n+1 to n jets (sometimes known as d_{n n+1}). If the number
126 /// of particles in the event is <= njets, the function returns 0.
127 double exclusive_dmerge (const int njets) const;
128
129 /// return the maximum of the dmin encountered during all recombinations
130 /// up to the one that led to an n-jet final state; identical to
131 /// exclusive_dmerge, except in cases where the dmin do not increase
132 /// monotonically.
133 double exclusive_dmerge_max (const int njets) const;
134
135 /// return the ymin corresponding to the recombination that went from
136 /// n+1 to n jets (sometimes known as y_{n n+1}).
137 double exclusive_ymerge (int njets) const {return exclusive_dmerge(njets) / Q2();}
138
139 /// same as exclusive_dmerge_max, but normalised to squared total energy
140 double exclusive_ymerge_max (int njets) const {return exclusive_dmerge_max(njets)/Q2();}
141
142 /// the number of exclusive jets at the given ycut
143 int n_exclusive_jets_ycut (double ycut) const {return n_exclusive_jets(ycut*Q2());}
144
145 /// the exclusive jets obtained at the given ycut
146 std::vector<PseudoJet> exclusive_jets_ycut (double ycut) const {
147 int njets = n_exclusive_jets_ycut(ycut);
148 return exclusive_jets(njets);
149 }
150
151
152 //int n_exclusive_jets (const PseudoJet & jet, const double dcut) const;
153
154 /// return a vector of all subjets of the current jet (in the sense
155 /// of the exclusive algorithm) that would be obtained when running
156 /// the algorithm with the given dcut.
157 ///
158 /// Time taken is O(m ln m), where m is the number of subjets that
159 /// are found. If m gets to be of order of the total number of
160 /// constituents in the jet, this could be substantially slower than
161 /// just getting that list of constituents.
162 std::vector<PseudoJet> exclusive_subjets (const PseudoJet & jet,
163 const double dcut) const;
164
165 /// return the size of exclusive_subjets(...); still n ln n with same
166 /// coefficient, but marginally more efficient than manually taking
167 /// exclusive_subjets.size()
168 int n_exclusive_subjets(const PseudoJet & jet,
169 const double dcut) const;
170
171 /// return the list of subjets obtained by unclustering the supplied
172 /// jet down to nsub subjets. Throws an error if there are fewer than
173 /// nsub particles in the jet.
174 ///
175 /// This requires nsub ln nsub time
176 std::vector<PseudoJet> exclusive_subjets (const PseudoJet & jet,
177 int nsub) const;
178
179 /// return the list of subjets obtained by unclustering the supplied
180 /// jet down to nsub subjets (or all constituents if there are fewer
181 /// than nsub).
182 ///
183 /// This requires nsub ln nsub time
184 std::vector<PseudoJet> exclusive_subjets_up_to (const PseudoJet & jet,
185 int nsub) const;
186
187 /// returns the dij that was present in the merging nsub+1 -> nsub
188 /// subjets inside this jet.
189 ///
190 /// Returns 0 if there were nsub or fewer constituents in the jet.
191 double exclusive_subdmerge(const PseudoJet & jet, int nsub) const;
192
193 /// returns the maximum dij that occurred in the whole event at the
194 /// stage that the nsub+1 -> nsub merge of subjets occurred inside
195 /// this jet.
196 ///
197 /// Returns 0 if there were nsub or fewer constituents in the jet.
198 double exclusive_subdmerge_max(const PseudoJet & jet, int nsub) const;
199
200 //std::vector<PseudoJet> exclusive_jets (const PseudoJet & jet,
201 // const int njets) const;
202 //double exclusive_dmerge (const PseudoJet & jet, const int njets) const;
203
204 /// returns the sum of all energies in the event (relevant mainly for e+e-)
205 double Q() const {return _Qtot;}
206 /// return Q()^2
207 double Q2() const {return _Qtot*_Qtot;}
208
209 /// returns true iff the object is included in the jet.
210 ///
211 /// NB: this is only sensible if the object is already registered
212 /// within the cluster sequence, so you cannot use it with an input
213 /// particle to the CS (since the particle won't have the history
214 /// index set properly).
215 ///
216 /// For nice clustering structures it should run in O(ln(N)) time
217 /// but in worst cases (certain cone plugins) it can take O(n) time,
218 /// where n is the number of particles in the jet.
219 bool object_in_jet(const PseudoJet & object, const PseudoJet & jet) const;
220
221 /// if the jet has parents in the clustering, it returns true
222 /// and sets parent1 and parent2 equal to them.
223 ///
224 /// if it has no parents it returns false and sets parent1 and
225 /// parent2 to zero
226 bool has_parents(const PseudoJet & jet, PseudoJet & parent1,
227 PseudoJet & parent2) const;
228
229 /// if the jet has a child then return true and give the child jet
230 /// otherwise return false and set the child to zero
231 bool has_child(const PseudoJet & jet, PseudoJet & child) const;
232
233 /// Version of has_child that sets a pointer to the child if the child
234 /// exists;
235 bool has_child(const PseudoJet & jet, const PseudoJet * & childp) const;
236
237 /// if this jet has a child (and so a partner) return true
238 /// and give the partner, otherwise return false and set the
239 /// partner to zero
240 bool has_partner(const PseudoJet & jet, PseudoJet & partner) const;
241
242
243 /// return a vector of the particles that make up jet
244 std::vector<PseudoJet> constituents (const PseudoJet & jet) const;
245
246
247 /// output the supplied vector of jets in a format that can be read
248 /// by an appropriate root script; the format is:
249 /// jet-n jet-px jet-py jet-pz jet-E
250 /// particle-n particle-rap particle-phi particle-pt
251 /// particle-n particle-rap particle-phi particle-pt
252 /// ...
253 /// #END
254 /// ... [i.e. above repeated]
255 void print_jets_for_root(const std::vector<PseudoJet> & jets,
256 std::ostream & ostr = std::cout) const;
257
258 /// print jets for root to the file labelled filename, with an
259 /// optional comment at the beginning
260 void print_jets_for_root(const std::vector<PseudoJet> & jets,
261 const std::string & filename,
262 const std::string & comment = "") const;
263
264// Not yet. Perhaps in a future release.
265// /// print out all inclusive jets with pt > ptmin
266// virtual void print_jets (const double ptmin=0.0) const;
267
268 /// add on to subjet_vector the constituents of jet (for internal use mainly)
269 void add_constituents (const PseudoJet & jet,
270 std::vector<PseudoJet> & subjet_vector) const;
271
272 /// return the enum value of the strategy used to cluster the event
273 inline Strategy strategy_used () const {return _strategy;}
274
275 /// return the name of the strategy used to cluster the event
276 std::string strategy_string () const {return strategy_string(_strategy);}
277
278 /// return the name of the strategy associated with the enum strategy_in
279 std::string strategy_string (Strategy strategy_in) const;
280
281
282 /// return a reference to the jet definition
283 const JetDefinition & jet_def() const {return _jet_def;}
284
285 /// by calling this routine you tell the ClusterSequence to delete
286 /// itself when all the Pseudojets associated with it have gone out
287 /// of scope.
288 ///
289 /// At the time you call this, there must be at least one jet or
290 /// other object outside the CS that is associated with the CS
291 /// (e.g. the result of inclusive_jets()).
292 ///
293 /// NB: after having made this call, the user is still allowed to
294 /// delete the CS. Jets associated with it will then simply not be
295 /// able to access their substructure after that point.
296 void delete_self_when_unused();
297
298 /// return true if the object has been told to delete itself
299 /// when unused
300 bool will_delete_self_when_unused() const {return _deletes_self_when_unused;}
301
302//std::shared_ptr: #ifdef FASTJET_HAVE_THREAD_SAFETY
303//std::shared_ptr: /// signals that a jet will no longer use the current CS (internal use only)
304//std::shared_ptr: void release_pseudojet(PseudoJet &jet) const;
305//std::shared_ptr: #else
306 /// tell the ClusterSequence it's about to be self deleted (internal use only)
307 void signal_imminent_self_deletion() const;
308//std::shared_ptr: #endif
309
310 /// returns the scale associated with a jet as required for this
311 /// clustering algorithm (kt^2 for the kt-algorithm, 1 for the
312 /// Cambridge algorithm). Intended mainly for internal use and not
313 /// valid for plugin algorithms.
314 double jet_scale_for_algorithm(const PseudoJet & jet) const;
315
316 ///
317
318 //----- next follow functions designed specifically for plugins, which
319 // may only be called when plugin_activated() returns true
320
321 /// record the fact that there has been a recombination between
322 /// jets()[jet_i] and jets()[jet_k], with the specified dij, and
323 /// return the index (newjet_k) allocated to the new jet, whose
324 /// momentum is assumed to be the 4-vector sum of that of jet_i and
325 /// jet_j
326 void plugin_record_ij_recombination(int jet_i, int jet_j, double dij,
327 int & newjet_k) {
328 assert(plugin_activated());
329 _do_ij_recombination_step(jet_i, jet_j, dij, newjet_k);
330 }
331
332 /// as for the simpler variant of plugin_record_ij_recombination,
333 /// except that the new jet is attributed the momentum and
334 /// user_index of newjet
335 void plugin_record_ij_recombination(int jet_i, int jet_j, double dij,
336 const PseudoJet & newjet,
337 int & newjet_k);
338
339 /// record the fact that there has been a recombination between
340 /// jets()[jet_i] and the beam, with the specified diB; when looking
341 /// for inclusive jets, any iB recombination will returned to the user
342 /// as a jet.
343 void plugin_record_iB_recombination(int jet_i, double diB) {
344 assert(plugin_activated());
345 _do_iB_recombination_step(jet_i, diB);
346 }
347
348 /// @ingroup extra_info
349 /// \class Extras
350 /// base class to store extra information that plugins may provide
351 ///
352 /// a class intended to serve as a base in case a plugin needs to
353 /// associate extra information with a ClusterSequence (see
354 /// SISConePlugin.* for an example).
355 class Extras {
356 public:
357 virtual ~Extras() {}
358 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";}
359 };
360
361 /// the plugin can associate some extra information with the
362 /// ClusterSequence object by calling this function. The
363 /// ClusterSequence takes ownership of the pointer (and
364 /// responsibility for deleting it when the CS gets deleted).
365 inline void plugin_associate_extras(Extras * extras_in) {
366 _extras.reset(extras_in);
367 }
368
369 /// the plugin can associate some extra information with the
370 /// ClusterSequence object by calling this function
371 ///
372 /// As of FJ v3.1, this is deprecated, in line with the deprecation
373 /// of auto_ptr in C++11
374#ifndef SWIG
375#ifdef FASTJET_HAVE_AUTO_PTR_INTERFACE
376 FASTJET_DEPRECATED_MSG("Please use ClusterSequence::plugin_associate_extras(Extras * extras_in)) instead",
377 inline void plugin_associate_extras(std::auto_ptr<Extras> extras_in)){
378 _extras.reset(extras_in.release());
379 }
380#endif
381#endif //SWIG
382
383 /// returns true when the plugin is allowed to run the show.
384 inline bool plugin_activated() const {return _plugin_activated;}
385
386 /// returns a pointer to the extras object (may be null)
387 const Extras * extras() const {return _extras.get();}
388
389 /// allows a plugin to run a templated clustering (nearest-neighbour heuristic)
390 ///
391 /// This has N^2 behaviour on "good" distance, but a worst case behaviour
392 /// of N^3 (and many algs trigger the worst case behaviour)
393 ///
394 ///
395 /// For more details on how this works, see GenBriefJet below
396 template<class GBJ> void plugin_simple_N2_cluster () {
397 assert(plugin_activated());
398 _simple_N2_cluster<GBJ>();
399 }
400
401
402public:
403//DEP /// set the default (static) jet finder across all current and future
404//DEP /// ClusterSequence objects -- deprecated and obsolescent (i.e. may be
405//DEP /// suppressed in a future release).
406//DEP static void set_jet_algorithm (JetAlgorithm jet_algorithm) {_default_jet_algorithm = jet_algorithm;}
407//DEP /// same as above for backward compatibility
408//DEP static void set_jet_finder (JetAlgorithm jet_algorithm) {_default_jet_algorithm = jet_algorithm;}
409
410
411 /// \ingroup extra_info
412 /// \struct history_element
413 /// a single element in the clustering history
414 ///
415 /// (see vector _history below).
417 int parent1; /// index in _history where first parent of this jet
418 /// was created (InexistentParent if this jet is an
419 /// original particle)
420
421 int parent2; /// index in _history where second parent of this jet
422 /// was created (InexistentParent if this jet is an
423 /// original particle); BeamJet if this history entry
424 /// just labels the fact that the jet has recombined
425 /// with the beam)
426
427 int child; /// index in _history where the current jet is
428 /// recombined with another jet to form its child. It
429 /// is Invalid if this jet does not further
430 /// recombine.
431
432 int jetp_index; /// index in the _jets vector where we will find the
433 /// PseudoJet object corresponding to this jet
434 /// (i.e. the jet created at this entry of the
435 /// history). NB: if this element of the history
436 /// corresponds to a beam recombination, then
437 /// jetp_index=Invalid.
438
439 double dij; /// the distance corresponding to the recombination
440 /// at this stage of the clustering.
441
442 double max_dij_so_far; /// the largest recombination distance seen
443 /// so far in the clustering history.
444 };
445
446 enum JetType {Invalid=-3, InexistentParent = -2, BeamJet = -1};
447
448 /// allow the user to access the internally stored _jets() array,
449 /// which contains both the initial particles and the various
450 /// intermediate and final stages of recombination.
451 ///
452 /// The first n_particles() entries are the original particles,
453 /// in the order in which they were supplied to the ClusterSequence
454 /// constructor. It can be useful to access them for example when
455 /// examining whether a given input object is part of a specific
456 /// jet, via the objects_in_jet(...) member function (which only takes
457 /// PseudoJets that are registered in the ClusterSequence).
458 ///
459 /// One of the other (internal uses) is related to the fact
460 /// because we don't seem to be able to access protected elements of
461 /// the class for an object that is not "this" (at least in case where
462 /// "this" is of a slightly different kind from the object, both
463 /// derived from ClusterSequence).
464 const std::vector<PseudoJet> & jets() const;
465
466 /// allow the user to access the raw internal history.
467 ///
468 /// This is present (as for jets()) in part so that protected
469 /// derived classes can access this information about other
470 /// ClusterSequences.
471 ///
472 /// A user who wishes to follow the details of the ClusterSequence
473 /// can also make use of this information (and should consult the
474 /// history_element documentation for more information), but should
475 /// be aware that these internal structures may evolve in future
476 /// FastJet versions.
477 const std::vector<history_element> & history() const;
478
479 /// returns the number of particles that were provided to the
480 /// clustering algorithm (helps the user find their way around the
481 /// history and jets objects if they weren't paying attention
482 /// beforehand).
483 unsigned int n_particles() const;
484
485 /// returns a vector of size n_particles() which indicates, for
486 /// each of the initial particles (in the order in which they were
487 /// supplied), which of the supplied jets it belongs to; if it does
488 /// not belong to any of the supplied jets, the index is set to -1;
489 std::vector<int> particle_jet_indices(const std::vector<PseudoJet> &) const;
490
491 /// routine that returns an order in which to read the history
492 /// such that clusterings that lead to identical jet compositions
493 /// but different histories (because of degeneracies in the
494 /// clustering order) will have matching constituents for each
495 /// matching entry in the unique_history_order.
496 ///
497 /// The order has the property that an entry's parents will always
498 /// appear prior to that entry itself.
499 ///
500 /// Roughly speaking the order is such that we first provide all
501 /// steps that lead to the final jet containing particle 1; then we
502 /// have the steps that lead to reconstruction of the jet containing
503 /// the next-lowest-numbered unclustered particle, etc...
504 /// [see GPS CCN28-12 for more info -- of course a full explanation
505 /// here would be better...]
506 std::vector<int> unique_history_order() const;
507
508 /// return the set of particles that have not been clustered. For
509 /// kt and cam/aachen algorithms this should always be null, but for
510 /// cone type algorithms it can be non-null;
511 std::vector<PseudoJet> unclustered_particles() const;
512
513 /// Return the list of pseudojets in the ClusterSequence that do not
514 /// have children (and are not among the inclusive jets). They may
515 /// result from a clustering step or may be one of the pseudojets
516 /// returned by unclustered_particles().
517 std::vector<PseudoJet> childless_pseudojets() const;
518
519 /// returns true if the object (jet or particle) is contained by (ie
520 /// belongs to) this cluster sequence.
521 ///
522 /// Tests performed: if thejet's interface is this cluster sequence
523 /// and its cluster history index is in a consistent range.
524 bool contains(const PseudoJet & object) const;
525
526 /// transfer the sequence contained in other_seq into our own;
527 /// any plugin "extras" contained in the from_seq will be lost
528 /// from there.
529 ///
530 /// It also sets the ClusterSequence pointers of the PseudoJets in
531 /// the history to point to this ClusterSequence
532 ///
533 /// When specified, the second argument is an action that will be
534 /// applied on every jets in the resulting ClusterSequence
535 void transfer_from_sequence(const ClusterSequence & from_seq,
536 const FunctionOfPseudoJet<PseudoJet> * action_on_jets = 0);
537
538 /// retrieve a shared pointer to the wrapper to this ClusterSequence
539 ///
540 /// this may turn useful if you want to track when this
541 /// ClusterSequence goes out of scope
543 return _structure_shared_ptr;
544 }
545
546 /// the structure type associated with a jet belonging to a ClusterSequence
548
549 /// This is the function that is automatically called during
550 /// clustering to print the FastJet banner. Only the first call to
551 /// this function will result in the printout of the banner. Users
552 /// may wish to call this function themselves, during the
553 /// initialization phase of their program, in order to ensure that
554 /// the banner appears before other output. This call will not
555 /// affect 3rd-party banners, e.g. those from plugins.
556 static void print_banner();
557
558 /// \cond internal_doc
559 // [this line must be left as is to hide the doxygen comment]
560 /// A call to this function modifies the stream used to print
561 /// banners (by default cout). If a null pointer is passed, banner
562 /// printout is suppressed. This affects all banners, including
563 /// those from plugins.
564 ///
565 /// Please note that if you distribute 3rd party code
566 /// that links with FastJet, that 3rd party code must not
567 /// use this call turn off the printing of FastJet banners
568 /// by default. This requirement reflects the spirit of
569 /// clause 2c of the GNU Public License (v2), under which
570 /// FastJet and its plugins are distributed.
571 static void set_fastjet_banner_stream(std::ostream * ostr) {_fastjet_banner_ostr = ostr;}
572 // [this line must be left as is to hide the doxygen comment]
573 /// \endcond
574
575 /// returns a pointer to the stream to be used to print banners
576 /// (cout by default). This function is used by plugins to determine
577 /// where to direct their banners. Plugins should properly handle
578 /// the case where the pointer is null.
579 static std::ostream * fastjet_banner_stream() {return _fastjet_banner_ostr;}
580
581private:
582
583 /// \cond internal_doc
584 /// contains the actual stream to use for banners
585#ifdef FASTJET_HAVE_LIMITED_THREAD_SAFETY
586 static std::atomic<std::ostream*> _fastjet_banner_ostr;
587#else
588 static std::ostream * _fastjet_banner_ostr;
589#endif // FASTJET_HAVE_LIMITED_THREAD_SAFETY
590 /// \endcond
591
592protected:
593//DEP static JetAlgorithm _default_jet_algorithm;
594 JetDefinition _jet_def;
595
596 /// transfer the vector<L> of input jets into our own vector<PseudoJet>
597 /// _jets (with some reserved space for future growth).
598 template<class L> void _transfer_input_jets(
599 const std::vector<L> & pseudojets);
600
601 /// This is what is called to do all the initialisation and
602 /// then run the clustering (may be called by various constructors).
603 /// It assumes _jets contains the momenta to be clustered.
604 void _initialise_and_run (const JetDefinition & jet_def,
605 const bool & writeout_combinations);
606
607 //// this performs the initialisation, minus the option-decanting
608 //// stage; for low multiplicity, initialising a few things in the
609 //// constructor, calling the decant_options_partial() and then this
610 //// is faster than going through _initialise_and_run.
611 void _initialise_and_run_no_decant();
612
613//DEP /// This is an alternative routine for initialising and running the
614//DEP /// clustering, provided for legacy purposes. The jet finder is that
615//DEP /// specified in the static member _default_jet_algorithm.
616//DEP void _initialise_and_run (const double R,
617//DEP const Strategy & strategy,
618//DEP const bool & writeout_combinations);
619
620 /// fills in the various member variables with "decanted" options from
621 /// the jet_definition and writeout_combinations variables
622 void _decant_options(const JetDefinition & jet_def,
623 const bool & writeout_combinations);
624
625 /// assuming that the jet definition, writeout_combinations and
626 /// _structure_shared_ptr have been set (e.g. in an initialiser list
627 /// in the constructor), it handles the remaining decanting of
628 /// options.
629 void _decant_options_partial();
630
631 /// fill out the history (and jet cross refs) related to the initial
632 /// set of jets (assumed already to have been "transferred"),
633 /// without any clustering
634 void _fill_initial_history();
635
636 /// carry out the recombination between the jets numbered jet_i and
637 /// jet_j, at distance scale dij; return the index newjet_k of the
638 /// result of the recombination of i and j.
639 void _do_ij_recombination_step(const int jet_i, const int jet_j,
640 const double dij, int & newjet_k);
641
642 /// carry out an recombination step in which _jets[jet_i] merges with
643 /// the beam,
644 void _do_iB_recombination_step(const int jet_i, const double diB);
645
646 /// every time a jet is added internally during clustering, this
647 /// should be called to set the jet's structure shared ptr to point
648 /// to the CS (and the count of internally associated objects is
649 /// also updated). This should not be called outside construction of
650 /// a CS object.
651 void _set_structure_shared_ptr(PseudoJet & j);
652
653 /// make sure that the CS's internal tally of the use count matches
654 /// that of the _structure_shared_ptr
655 void _update_structure_use_count();
656
657 /// returns a suggestion for the best strategy to use on event
658 /// multiplicity, algorithm, R, etc.
659 Strategy _best_strategy() const;
660
661 /// \if internal_doc
662 /// \class _Parabola
663 /// returns c*(a*R**2 + b*R + 1);
664 /// Written as a class in case we want to give names to different
665 /// parabolas
666 /// \endif
667 class _Parabola {
668 public:
669 _Parabola(double a, double b, double c) : _a(a), _b(b), _c(c) {}
670 inline double operator()(const double R) const {return _c*(_a*R*R + _b*R + 1);}
671 private:
672 double _a, _b, _c;
673 };
674
675 /// \if internal_doc
676 /// \class _Line
677 /// operator()(R) returns a*R+b;
678 /// \endif
679 class _Line {
680 public:
681 _Line(double a, double b) : _a(a), _b(b) {}
682 inline double operator()(const double R) const {return _a*R + _b;}
683 private:
684 double _a, _b;
685 };
686
687 /// This contains the physical PseudoJets; for each PseudoJet one
688 /// can find the corresponding position in the _history by looking
689 /// at _jets[i].cluster_hist_index().
690 std::vector<PseudoJet> _jets;
691
692
693 /// this vector will contain the branching history; for each stage,
694 /// _history[i].jetp_index indicates where to look in the _jets
695 /// vector to get the physical PseudoJet.
696 std::vector<history_element> _history;
697
698 /// set subhist to be a set pointers to history entries corresponding to the
699 /// subjets of this jet; one stops going working down through the
700 /// subjets either when
701 /// - there is no further to go
702 /// - one has found maxjet entries
703 /// - max_dij_so_far <= dcut
704 /// By setting maxjet=0 one can use just dcut; by setting dcut<0
705 /// one can use jet maxjet
706 void get_subhist_set(std::set<const history_element*> & subhist,
707 const PseudoJet & jet, double dcut, int maxjet) const;
708
709 bool _writeout_combinations;
710 int _initial_n;
711 double _Rparam, _R2, _invR2;
712 double _Qtot;
713 Strategy _strategy;
714 JetAlgorithm _jet_algorithm;
715
716 SharedPtr<PseudoJetStructureBase> _structure_shared_ptr; //< will actually be of type ClusterSequenceStructure
717 int _structure_use_count_after_construction; //< info of use when CS handles its own memory
718 /// if true then the CS will delete itself when the last external
719 /// object referring to it disappears. It is mutable so as to ensure
720 /// that signal_imminent_self_deletion() [const] can make relevant
721 /// changes.
722#ifdef FASTJET_HAVE_THREAD_SAFETY
723 mutable std::atomic<bool> _deletes_self_when_unused;
724#else
726#endif
727
728 private:
729
730 bool _plugin_activated;
731 SharedPtr<Extras> _extras; // things the plugin might want to add
732
733 void _really_dumb_cluster ();
734 void _delaunay_cluster ();
735 //void _simple_N2_cluster ();
736 template<class BJ> void _simple_N2_cluster ();
737 void _tiled_N2_cluster ();
738 void _faster_tiled_N2_cluster ();
739
740 //
741 void _minheap_faster_tiled_N2_cluster();
742
743 // things needed specifically for Cambridge with Chan's 2D closest
744 // pairs method
745 void _CP2DChan_cluster();
746 void _CP2DChan_cluster_2pi2R ();
747 void _CP2DChan_cluster_2piMultD ();
748 void _CP2DChan_limited_cluster(double D);
749 void _do_Cambridge_inclusive_jets();
750
751 // NSqrtN method for C/A
752 void _fast_NsqrtN_cluster();
753
754 void _add_step_to_history( //const int step_number,
755 const int parent1,
756 const int parent2, const int jetp_index,
757 const double dij);
758
759 /// internal routine associated with the construction of the unique
760 /// history order (following children in the tree)
761 void _extract_tree_children(int pos, std::valarray<bool> &,
762 const std::valarray<int> &, std::vector<int> &) const;
763
764 /// internal routine associated with the construction of the unique
765 /// history order (following parents in the tree)
766 void _extract_tree_parents (int pos, std::valarray<bool> &,
767 const std::valarray<int> &, std::vector<int> &) const;
768
769
770 // these will be useful shorthands in the Voronoi-based code
771 typedef std::pair<int,int> TwoVertices;
772 typedef std::pair<double,TwoVertices> DijEntry;
773 typedef std::multimap<double,TwoVertices> DistMap;
774
775 /// currently used only in the Voronoi based code
776 void _add_ktdistance_to_map(const int ii,
777 DistMap & DijMap,
778 const DynamicNearestNeighbours * DNN);
779
780
781 /// will be set by default to be true for the first run
782 static thread_safety_helpers::FirstTimeTrue _first_time;
783
784 /// manage warnings related to exclusive jets access
785 static LimitedWarning _exclusive_warnings;
786
787 /// the limited warning member for notification of user that
788 /// their requested strategy has been overridden (usually because
789 /// they have R>2pi and not all strategies work then)
790 static LimitedWarning _changed_strategy_warning;
791
792 //----------------------------------------------------------------------
793 /// the fundamental structure which contains the minimal info about
794 /// a jet, as needed for our plain N^2 algorithm -- the idea is to
795 /// put all info that will be accessed N^2 times into an array of
796 /// BriefJets...
797 struct BriefJet {
798 double eta, phi, kt2, NN_dist;
799 BriefJet * NN;
800 int _jets_index;
801 };
802
803 /// structure analogous to BriefJet, but with the extra information
804 /// needed for dealing with tiles
805 class TiledJet {
806 public:
807 double eta, phi, kt2, NN_dist;
808 TiledJet * NN, *previous, * next;
809 int _jets_index, tile_index, diJ_posn;
810 // routines that are useful in the minheap version of tiled
811 // clustering ("misuse" the otherwise unused diJ_posn, so as
812 // to indicate whether jets need to have their minheap entries
813 // updated).
814 inline void label_minheap_update_needed() {diJ_posn = 1;}
815 inline void label_minheap_update_done() {diJ_posn = 0;}
816 inline bool minheap_update_needed() const {return diJ_posn==1;}
817 };
818
819 //-- some of the functions that follow are templates and will work
820 //as well for briefjet and tiled jets
821
822 /// set the kinematic and labelling info for jeta so that it corresponds
823 /// to _jets[_jets_index]
824 template <class J> void _bj_set_jetinfo( J * const jet,
825 const int _jets_index) const;
826
827 /// "remove" this jet, which implies updating links of neighbours and
828 /// perhaps modifying the tile structure
829 void _bj_remove_from_tiles( TiledJet * const jet) const;
830
831 /// return the distance between two BriefJet objects
832 template <class J> double _bj_dist(const J * const jeta,
833 const J * const jetb) const;
834
835 // return the diJ (multiplied by _R2) for this jet assuming its NN
836 // info is correct
837 template <class J> double _bj_diJ(const J * const jeta) const;
838
839 /// for testing purposes only: if in the range head--tail-1 there is a
840 /// a jet which corresponds to hist_index in the history, then
841 /// return a pointer to that jet; otherwise return tail.
842 template <class J> inline J * _bj_of_hindex(
843 const int hist_index,
844 J * const head, J * const tail)
845 const {
846 J * res;
847 for(res = head; res<tail; res++) {
848 if (_jets[res->_jets_index].cluster_hist_index() == hist_index) {break;}
849 }
850 return res;
851 }
852
853
854 //-- remaining functions are different in various cases, so we
855 // will use templates but are not sure if they're useful...
856
857 /// updates (only towards smaller distances) the NN for jeta without checking
858 /// whether in the process jeta itself might be a new NN of one of
859 /// the jets being scanned -- span the range head to tail-1 with
860 /// assumption that jeta is not contained in that range
861 template <class J> void _bj_set_NN_nocross(J * const jeta,
862 J * const head, const J * const tail) const;
863
864 /// reset the NN for jeta and DO check whether in the process jeta
865 /// itself might be a new NN of one of the jets being scanned --
866 /// span the range head to tail-1 with assumption that jeta is not
867 /// contained in that range
868 template <class J> void _bj_set_NN_crosscheck(J * const jeta,
869 J * const head, const J * const tail) const;
870
871
872
873 /// number of neighbours that a tile will have (rectangular geometry
874 /// gives 9 neighbours).
875 static const int n_tile_neighbours = 9;
876 //----------------------------------------------------------------------
877 /// The fundamental structures to be used for the tiled N^2 algorithm
878 /// (see CCN27-44 for some discussion of pattern of tiling)
879 struct Tile {
880 /// pointers to neighbouring tiles, including self
881 Tile * begin_tiles[n_tile_neighbours];
882 /// neighbouring tiles, excluding self
883 Tile ** surrounding_tiles;
884 /// half of neighbouring tiles, no self
885 Tile ** RH_tiles;
886 /// just beyond end of tiles
887 Tile ** end_tiles;
888 /// start of list of BriefJets contained in this tile
889 TiledJet * head;
890 /// sometimes useful to be able to tag a tile
891 bool tagged;
892 };
893 std::vector<Tile> _tiles;
894 double _tiles_eta_min, _tiles_eta_max;
895 double _tile_size_eta, _tile_size_phi;
896 int _n_tiles_phi,_tiles_ieta_min,_tiles_ieta_max;
897
898 // reasonably robust return of tile index given ieta and iphi, in particular
899 // it works even if iphi is negative
900 inline int _tile_index (int ieta, int iphi) const {
901 // note that (-1)%n = -1 so that we have to add _n_tiles_phi
902 // before performing modulo operation
903 return (ieta-_tiles_ieta_min)*_n_tiles_phi
904 + (iphi+_n_tiles_phi) % _n_tiles_phi;
905 }
906
907 // routines for tiled case, including some overloads of the plain
908 // BriefJet cases
909 int _tile_index(const double eta, const double phi) const;
910 void _tj_set_jetinfo ( TiledJet * const jet, const int _jets_index);
911 void _bj_remove_from_tiles(TiledJet * const jet);
912 void _initialise_tiles();
913 void _print_tiles(TiledJet * briefjets ) const;
914 void _add_neighbours_to_tile_union(const int tile_index,
915 std::vector<int> & tile_union, int & n_near_tiles) const;
916 void _add_untagged_neighbours_to_tile_union(const int tile_index,
917 std::vector<int> & tile_union, int & n_near_tiles);
918
919 //----------------------------------------------------------------------
920 /// fundamental structure for e+e- clustering
921 struct EEBriefJet {
922 double NN_dist; // obligatorily present
923 double kt2; // obligatorily present == E^2 in general
924 EEBriefJet * NN; // must be present too
925 int _jets_index; // must also be present!
926 //...........................................................
927 double nx, ny, nz; // our internal storage for fast distance calcs
928 };
929
930 /// to help instantiation (fj 2.4.0; did not quite work on gcc 33 and os x 10.3?)
931 //void _dummy_N2_cluster_instantiation();
932
933
934 /// to avoid issues with template instantiation (OS X 10.3, gcc 3.3)
935 void _simple_N2_cluster_BriefJet();
936 /// to avoid issues with template instantiation (OS X 10.3, gcc 3.3)
937 void _simple_N2_cluster_EEBriefJet();
938};
939
940
941//**********************************************************************
942//************** START OF INLINE MATERIAL ******************
943//**********************************************************************
944
945
946//----------------------------------------------------------------------
947// Transfer the initial jets into our internal structure
948template<class L> void ClusterSequence::_transfer_input_jets(
949 const std::vector<L> & pseudojets) {
950
951 // this will ensure that we can point to jets without difficulties
952 // arising.
953 _jets.reserve(pseudojets.size()*2);
954
955 // insert initial jets this way so that any type L that can be
956 // converted to a pseudojet will work fine (basically PseudoJet
957 // and any type that has [] subscript access to the momentum
958 // components, such as CLHEP HepLorentzVector).
959 for (unsigned int i = 0; i < pseudojets.size(); i++) {
960 _jets.push_back(pseudojets[i]);}
961
962}
963
964// //----------------------------------------------------------------------
965// // initialise from some generic type... Has to be made available
966// // here in order for it the template aspect of it to work...
967// template<class L> ClusterSequence::ClusterSequence (
968// const std::vector<L> & pseudojets,
969// const double R,
970// const Strategy & strategy,
971// const bool & writeout_combinations) {
972//
973// // transfer the initial jets (type L) into our own array
974// _transfer_input_jets(pseudojets);
975//
976// // run the clustering
977// _initialise_and_run(R,strategy,writeout_combinations);
978// }
979
980
981//----------------------------------------------------------------------
982/// constructor of a jet-clustering sequence from a vector of
983/// four-momenta, with the jet definition specified by jet_def
984template<class L> ClusterSequence::ClusterSequence (
985 const std::vector<L> & pseudojets,
986 const JetDefinition & jet_def_in,
987 const bool & writeout_combinations) :
988 _jet_def(jet_def_in), _writeout_combinations(writeout_combinations),
989 _structure_shared_ptr(new ClusterSequenceStructure(this))
990{
991
992 // transfer the initial jets (type L) into our own array
993 _transfer_input_jets(pseudojets);
994
995 // transfer the remaining options
997
998 // run the clustering
999 _initialise_and_run_no_decant();
1000}
1001
1002
1003inline const std::vector<PseudoJet> & ClusterSequence::jets () const {
1004 return _jets;
1005}
1006
1007inline const std::vector<ClusterSequence::history_element> & ClusterSequence::history () const {
1008 return _history;
1009}
1010
1011inline unsigned int ClusterSequence::n_particles() const {return _initial_n;}
1012
1013//----------------------------------------------------------------------
1014// implementation of JetDefinition::operator() is here to avoid nasty
1015// issues of order of implementations and includes
1016#ifndef __CINT__
1017template<class L>
1018std::vector<PseudoJet> JetDefinition::operator()(const std::vector<L> & particles) const {
1019 // create a new cluster sequence
1020 ClusterSequence * cs = new ClusterSequence(particles, *this);
1021
1022 // get the jets, and sort them according to whether the algorithm
1023 // is spherical or not
1024 std::vector<PseudoJet> jets;
1025 if (is_spherical()) {
1026 jets = sorted_by_E(cs->inclusive_jets());
1027 } else {
1028 jets = sorted_by_pt(cs->inclusive_jets());
1029 }
1030
1031 // make sure the ClusterSequence gets deleted once it's no longer
1032 // needed
1033 if (jets.size() != 0) {
1035 } else {
1036 delete cs;
1037 }
1038
1039 return jets;
1040}
1041#endif // __CINT__
1042
1043
1044//----------------------------------------------------------------------
1045template <class J> inline void ClusterSequence::_bj_set_jetinfo(
1046 J * const jetA, const int _jets_index) const {
1047 jetA->eta = _jets[_jets_index].rap();
1048 jetA->phi = _jets[_jets_index].phi_02pi();
1049 jetA->kt2 = jet_scale_for_algorithm(_jets[_jets_index]);
1050 jetA->_jets_index = _jets_index;
1051 // initialise NN info as well
1052 jetA->NN_dist = _R2;
1053 jetA->NN = NULL;
1054}
1055
1056
1057
1058
1059//----------------------------------------------------------------------
1060template <class J> inline double ClusterSequence::_bj_dist(
1061 const J * const jetA, const J * const jetB) const {
1062 //#define FASTJET_NEW_DELTA_PHI
1063#ifndef FASTJET_NEW_DELTA_PHI
1064 //GPS+MC old version of Delta phi calculation
1065 double dphi = std::abs(jetA->phi - jetB->phi);
1066 double deta = (jetA->eta - jetB->eta);
1067 if (dphi > pi) {dphi = twopi - dphi;}
1068#else
1069 //GPS+MC testing for 2015-02-faster-deltaR2
1070 double dphi = pi-std::abs(pi-std::abs(jetA->phi - jetB->phi));
1071 double deta = (jetA->eta - jetB->eta);
1072#endif
1073 return dphi*dphi + deta*deta;
1074}
1075
1076//----------------------------------------------------------------------
1077template <class J> inline double ClusterSequence::_bj_diJ(const J * const jet) const {
1078 double kt2 = jet->kt2;
1079 if (jet->NN != NULL) {if (jet->NN->kt2 < kt2) {kt2 = jet->NN->kt2;}}
1080 return jet->NN_dist * kt2;
1081}
1082
1083
1084//----------------------------------------------------------------------
1085// set the NN for jet without checking whether in the process you might
1086// have discovered a new nearest neighbour for another jet
1087template <class J> inline void ClusterSequence::_bj_set_NN_nocross(
1088 J * const jet, J * const head, const J * const tail) const {
1089 double NN_dist = _R2;
1090 J * NN = NULL;
1091 if (head < jet) {
1092 for (J * jetB = head; jetB != jet; jetB++) {
1093 double dist = _bj_dist(jet,jetB);
1094 if (dist < NN_dist) {
1095 NN_dist = dist;
1096 NN = jetB;
1097 }
1098 }
1099 }
1100 if (tail > jet) {
1101 for (J * jetB = jet+1; jetB != tail; jetB++) {
1102 double dist = _bj_dist(jet,jetB);
1103 if (dist < NN_dist) {
1104 NN_dist = dist;
1105 NN = jetB;
1106 }
1107 }
1108 }
1109 jet->NN = NN;
1110 jet->NN_dist = NN_dist;
1111}
1112
1113
1114//----------------------------------------------------------------------
1115template <class J> inline void ClusterSequence::_bj_set_NN_crosscheck(J * const jet,
1116 J * const head, const J * const tail) const {
1117 double NN_dist = _R2;
1118 J * NN = NULL;
1119 for (J * jetB = head; jetB != tail; jetB++) {
1120 double dist = _bj_dist(jet,jetB);
1121 if (dist < NN_dist) {
1122 NN_dist = dist;
1123 NN = jetB;
1124 }
1125 if (dist < jetB->NN_dist) {
1126 jetB->NN_dist = dist;
1127 jetB->NN = jet;
1128 }
1129 }
1130 jet->NN = NN;
1131 jet->NN_dist = NN_dist;
1132}
1133
1134FASTJET_END_NAMESPACE
1135
1136#endif // __FASTJET_CLUSTERSEQUENCE_HH__
Contains any information related to the clustering that should be directly accessible to PseudoJet.
base class to store extra information that plugins may provide
deals with clustering
bool _deletes_self_when_unused
if true then the CS will delete itself when the last external object referring to it disappears.
const SharedPtr< PseudoJetStructureBase > & structure_shared_ptr() const
retrieve a shared pointer to the wrapper to this ClusterSequence
bool plugin_activated() const
returns true when the plugin is allowed to run the show.
std::vector< history_element > _history
this vector will contain the branching history; for each stage, _history[i].jetp_index indicates wher...
ClusterSequenceStructure StructureType
the structure type associated with a jet belonging to a ClusterSequence
void plugin_simple_N2_cluster()
allows a plugin to run a templated clustering (nearest-neighbour heuristic)
const JetDefinition & jet_def() const
return a reference to the jet definition
int n_exclusive_jets_ycut(double ycut) const
the number of exclusive jets at the given ycut
const std::vector< history_element > & history() const
allow the user to access the raw internal history.
void plugin_associate_extras(std::auto_ptr< Extras > extras_in)
the plugin can associate some extra information with the ClusterSequence object by calling this funct...
void plugin_record_iB_recombination(int jet_i, double diB)
record the fact that there has been a recombination between jets()[jet_i] and the beam,...
void delete_self_when_unused()
by calling this routine you tell the ClusterSequence to delete itself when all the Pseudojets associa...
unsigned int n_particles() const
returns the number of particles that were provided to the clustering algorithm (helps the user find t...
std::vector< PseudoJet > inclusive_jets(const double ptmin=0.0) const
return a vector of all jets (in the sense of the inclusive algorithm) with pt >= ptmin.
std::string strategy_string() const
return the name of the strategy used to cluster the event
void _decant_options_partial()
assuming that the jet definition, writeout_combinations and _structure_shared_ptr have been set (e....
std::vector< PseudoJet > _jets
This contains the physical PseudoJets; for each PseudoJet one can find the corresponding position in ...
bool will_delete_self_when_unused() const
return true if the object has been told to delete itself when unused
void print_jets_for_root(const std::vector< PseudoJet > &jets, std::ostream &ostr=std::cout) const
output the supplied vector of jets in a format that can be read by an appropriate root script; the fo...
const std::vector< PseudoJet > & jets() const
allow the user to access the internally stored _jets() array, which contains both the initial particl...
ClusterSequence(const ClusterSequence &cs)
copy constructor for a ClusterSequence
double exclusive_ymerge(int njets) const
return the ymin corresponding to the recombination that went from n+1 to n jets (sometimes known as y...
double exclusive_ymerge_max(int njets) const
same as exclusive_dmerge_max, but normalised to squared total energy
void plugin_associate_extras(Extras *extras_in)
the plugin can associate some extra information with the ClusterSequence object by calling this funct...
double Q() const
returns the sum of all energies in the event (relevant mainly for e+e-)
ClusterSequence()
default constructor
std::vector< PseudoJet > exclusive_jets_ycut(double ycut) const
the exclusive jets obtained at the given ycut
double Q2() const
return Q()^2
static std::ostream * fastjet_banner_stream()
returns a pointer to the stream to be used to print banners (cout by default).
void _transfer_input_jets(const std::vector< L > &pseudojets)
transfer the vector<L> of input jets into our own vector<PseudoJet> _jets (with some reserved space f...
double jet_scale_for_algorithm(const PseudoJet &jet) const
returns the scale associated with a jet as required for this clustering algorithm (kt^2 for the kt-al...
const Extras * extras() const
returns a pointer to the extras object (may be null)
void plugin_record_ij_recombination(int jet_i, int jet_j, double dij, int &newjet_k)
record the fact that there has been a recombination between jets()[jet_i] and jets()[jet_k],...
Strategy strategy_used() const
return the enum value of the strategy used to cluster the event
class that is intended to hold a full definition of the jet clusterer
bool is_spherical() const
returns true if the jet definition involves an algorithm intended for use on a spherical geometry (e....
std::vector< PseudoJet > operator()(const std::vector< L > &particles) const
cluster the supplied particles and returns a vector of resulting jets, sorted by pt (or energy in the...
class to provide facilities for giving warnings up to some maximum number of times and to provide glo...
Class to contain pseudojets, including minimal information of use to jet-clustering routines.
Definition: PseudoJet.hh:68
An implementation of shared pointers that is broadly similar to C++11 shared_ptr (https://en....
Definition: SharedPtr.hh:341
structure analogous to BriefJet, but with the extra information needed for dealing with tiles
provides an object wich will return "true" the first time () is called and false afterwards
Strategy
the various options for the algorithmic strategy to adopt in clustering events with kt and cambridge ...
vector< PseudoJet > sorted_by_E(const vector< PseudoJet > &jets)
return a vector of jets sorted into decreasing energy
Definition: PseudoJet.cc:898
JetAlgorithm
the various families of jet-clustering algorithm
vector< PseudoJet > sorted_by_pt(const vector< PseudoJet > &jets)
return a vector of jets sorted into decreasing kt2
Definition: PseudoJet.cc:882
a single element in the clustering history
int jetp_index
index in _history where the current jet is recombined with another jet to form its child.
int parent2
index in _history where first parent of this jet was created (InexistentParent if this jet is an orig...
int child
index in _history where second parent of this jet was created (InexistentParent if this jet is an ori...
double max_dij_so_far
the distance corresponding to the recombination at this stage of the clustering.
double dij
index in the _jets vector where we will find the