FastJet  3.3.3
ClusterSequence.hh
1 #ifndef __FASTJET_CLUSTERSEQUENCE_HH__
2 #define __FASTJET_CLUSTERSEQUENCE_HH__
3 
4 //FJSTARTHEADER
5 // $Id: ClusterSequence.hh 4420 2019-11-29 09:28:20Z soyez $
6 //
7 // Copyright (c) 2005-2019, 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 
51 #include "fastjet/internal/deprecated.hh"
52 
53 FASTJET_BEGIN_NAMESPACE // defined in fastjet/internal/base.hh
54 
55 
56 // forward declarations
57 class ClusterSequenceStructure;
58 class 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  /// tell the ClusterSequence it's about to be self deleted (internal use only)
303  void signal_imminent_self_deletion() const;
304 
305  /// returns the scale associated with a jet as required for this
306  /// clustering algorithm (kt^2 for the kt-algorithm, 1 for the
307  /// Cambridge algorithm). Intended mainly for internal use and not
308  /// valid for plugin algorithms.
309  double jet_scale_for_algorithm(const PseudoJet & jet) const;
310 
311  ///
312 
313  //----- next follow functions designed specifically for plugins, which
314  // may only be called when plugin_activated() returns true
315 
316  /// record the fact that there has been a recombination between
317  /// jets()[jet_i] and jets()[jet_k], with the specified dij, and
318  /// return the index (newjet_k) allocated to the new jet, whose
319  /// momentum is assumed to be the 4-vector sum of that of jet_i and
320  /// jet_j
321  void plugin_record_ij_recombination(int jet_i, int jet_j, double dij,
322  int & newjet_k) {
323  assert(plugin_activated());
324  _do_ij_recombination_step(jet_i, jet_j, dij, newjet_k);
325  }
326 
327  /// as for the simpler variant of plugin_record_ij_recombination,
328  /// except that the new jet is attributed the momentum and
329  /// user_index of newjet
330  void plugin_record_ij_recombination(int jet_i, int jet_j, double dij,
331  const PseudoJet & newjet,
332  int & newjet_k);
333 
334  /// record the fact that there has been a recombination between
335  /// jets()[jet_i] and the beam, with the specified diB; when looking
336  /// for inclusive jets, any iB recombination will returned to the user
337  /// as a jet.
338  void plugin_record_iB_recombination(int jet_i, double diB) {
339  assert(plugin_activated());
340  _do_iB_recombination_step(jet_i, diB);
341  }
342 
343  /// @ingroup extra_info
344  /// \class Extras
345  /// base class to store extra information that plugins may provide
346  ///
347  /// a class intended to serve as a base in case a plugin needs to
348  /// associate extra information with a ClusterSequence (see
349  /// SISConePlugin.* for an example).
350  class Extras {
351  public:
352  virtual ~Extras() {}
353  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";}
354  };
355 
356  /// the plugin can associate some extra information with the
357  /// ClusterSequence object by calling this function. The
358  /// ClusterSequence takes ownership of the pointer (and
359  /// responsibility for deleting it when the CS gets deleted).
360  inline void plugin_associate_extras(Extras * extras_in) {
361  _extras.reset(extras_in);
362  }
363 
364  /// the plugin can associate some extra information with the
365  /// ClusterSequence object by calling this function
366  ///
367  /// As of FJ v3.1, this is deprecated, in line with the deprecation
368  /// of auto_ptr in C++11
369 #ifndef SWIG
370 #ifdef FASTJET_HAVE_AUTO_PTR_INTERFACE
371  FASTJET_DEPRECATED_MSG("Please use ClusterSequence::plugin_associate_extras(Extras * extras_in)) instead")
372  inline void plugin_associate_extras(std::auto_ptr<Extras> extras_in){
373  _extras.reset(extras_in.release());
374  }
375 #endif
376 #endif //SWIG
377 
378  /// returns true when the plugin is allowed to run the show.
379  inline bool plugin_activated() const {return _plugin_activated;}
380 
381  /// returns a pointer to the extras object (may be null)
382  const Extras * extras() const {return _extras.get();}
383 
384  /// allows a plugin to run a templated clustering (nearest-neighbour heuristic)
385  ///
386  /// This has N^2 behaviour on "good" distance, but a worst case behaviour
387  /// of N^3 (and many algs trigger the worst case behaviour)
388  ///
389  ///
390  /// For more details on how this works, see GenBriefJet below
391  template<class GBJ> void plugin_simple_N2_cluster () {
392  assert(plugin_activated());
393  _simple_N2_cluster<GBJ>();
394  }
395 
396 
397 public:
398 //DEP /// set the default (static) jet finder across all current and future
399 //DEP /// ClusterSequence objects -- deprecated and obsolescent (i.e. may be
400 //DEP /// suppressed in a future release).
401 //DEP static void set_jet_algorithm (JetAlgorithm jet_algorithm) {_default_jet_algorithm = jet_algorithm;}
402 //DEP /// same as above for backward compatibility
403 //DEP static void set_jet_finder (JetAlgorithm jet_algorithm) {_default_jet_algorithm = jet_algorithm;}
404 
405 
406  /// \ingroup extra_info
407  /// \struct history_element
408  /// a single element in the clustering history
409  ///
410  /// (see vector _history below).
412  int parent1; /// index in _history where first parent of this jet
413  /// was created (InexistentParent if this jet is an
414  /// original particle)
415 
416  int parent2; /// index in _history where second parent of this jet
417  /// was created (InexistentParent if this jet is an
418  /// original particle); BeamJet if this history entry
419  /// just labels the fact that the jet has recombined
420  /// with the beam)
421 
422  int child; /// index in _history where the current jet is
423  /// recombined with another jet to form its child. It
424  /// is Invalid if this jet does not further
425  /// recombine.
426 
427  int jetp_index; /// index in the _jets vector where we will find the
428  /// PseudoJet object corresponding to this jet
429  /// (i.e. the jet created at this entry of the
430  /// history). NB: if this element of the history
431  /// corresponds to a beam recombination, then
432  /// jetp_index=Invalid.
433 
434  double dij; /// the distance corresponding to the recombination
435  /// at this stage of the clustering.
436 
437  double max_dij_so_far; /// the largest recombination distance seen
438  /// so far in the clustering history.
439  };
440 
441  enum JetType {Invalid=-3, InexistentParent = -2, BeamJet = -1};
442 
443  /// allow the user to access the internally stored _jets() array,
444  /// which contains both the initial particles and the various
445  /// intermediate and final stages of recombination.
446  ///
447  /// The first n_particles() entries are the original particles,
448  /// in the order in which they were supplied to the ClusterSequence
449  /// constructor. It can be useful to access them for example when
450  /// examining whether a given input object is part of a specific
451  /// jet, via the objects_in_jet(...) member function (which only takes
452  /// PseudoJets that are registered in the ClusterSequence).
453  ///
454  /// One of the other (internal uses) is related to the fact
455  /// because we don't seem to be able to access protected elements of
456  /// the class for an object that is not "this" (at least in case where
457  /// "this" is of a slightly different kind from the object, both
458  /// derived from ClusterSequence).
459  const std::vector<PseudoJet> & jets() const;
460 
461  /// allow the user to access the raw internal history.
462  ///
463  /// This is present (as for jets()) in part so that protected
464  /// derived classes can access this information about other
465  /// ClusterSequences.
466  ///
467  /// A user who wishes to follow the details of the ClusterSequence
468  /// can also make use of this information (and should consult the
469  /// history_element documentation for more information), but should
470  /// be aware that these internal structures may evolve in future
471  /// FastJet versions.
472  const std::vector<history_element> & history() const;
473 
474  /// returns the number of particles that were provided to the
475  /// clustering algorithm (helps the user find their way around the
476  /// history and jets objects if they weren't paying attention
477  /// beforehand).
478  unsigned int n_particles() const;
479 
480  /// returns a vector of size n_particles() which indicates, for
481  /// each of the initial particles (in the order in which they were
482  /// supplied), which of the supplied jets it belongs to; if it does
483  /// not belong to any of the supplied jets, the index is set to -1;
484  std::vector<int> particle_jet_indices(const std::vector<PseudoJet> &) const;
485 
486  /// routine that returns an order in which to read the history
487  /// such that clusterings that lead to identical jet compositions
488  /// but different histories (because of degeneracies in the
489  /// clustering order) will have matching constituents for each
490  /// matching entry in the unique_history_order.
491  ///
492  /// The order has the property that an entry's parents will always
493  /// appear prior to that entry itself.
494  ///
495  /// Roughly speaking the order is such that we first provide all
496  /// steps that lead to the final jet containing particle 1; then we
497  /// have the steps that lead to reconstruction of the jet containing
498  /// the next-lowest-numbered unclustered particle, etc...
499  /// [see GPS CCN28-12 for more info -- of course a full explanation
500  /// here would be better...]
501  std::vector<int> unique_history_order() const;
502 
503  /// return the set of particles that have not been clustered. For
504  /// kt and cam/aachen algorithms this should always be null, but for
505  /// cone type algorithms it can be non-null;
506  std::vector<PseudoJet> unclustered_particles() const;
507 
508  /// Return the list of pseudojets in the ClusterSequence that do not
509  /// have children (and are not among the inclusive jets). They may
510  /// result from a clustering step or may be one of the pseudojets
511  /// returned by unclustered_particles().
512  std::vector<PseudoJet> childless_pseudojets() const;
513 
514  /// returns true if the object (jet or particle) is contained by (ie
515  /// belongs to) this cluster sequence.
516  ///
517  /// Tests performed: if thejet's interface is this cluster sequence
518  /// and its cluster history index is in a consistent range.
519  bool contains(const PseudoJet & object) const;
520 
521  /// transfer the sequence contained in other_seq into our own;
522  /// any plugin "extras" contained in the from_seq will be lost
523  /// from there.
524  ///
525  /// It also sets the ClusterSequence pointers of the PseudoJets in
526  /// the history to point to this ClusterSequence
527  ///
528  /// When specified, the second argument is an action that will be
529  /// applied on every jets in the resulting ClusterSequence
530  void transfer_from_sequence(const ClusterSequence & from_seq,
531  const FunctionOfPseudoJet<PseudoJet> * action_on_jets = 0);
532 
533  /// retrieve a shared pointer to the wrapper to this ClusterSequence
534  ///
535  /// this may turn useful if you want to track when this
536  /// ClusterSequence goes out of scope
538  return _structure_shared_ptr;
539  }
540 
541  /// the structure type associated with a jet belonging to a ClusterSequence
543 
544  /// This is the function that is automatically called during
545  /// clustering to print the FastJet banner. Only the first call to
546  /// this function will result in the printout of the banner. Users
547  /// may wish to call this function themselves, during the
548  /// initialization phase of their program, in order to ensure that
549  /// the banner appears before other output. This call will not
550  /// affect 3rd-party banners, e.g. those from plugins.
551  static void print_banner();
552 
553  /// \cond internal_doc
554  // [this line must be left as is to hide the doxygen comment]
555  /// A call to this function modifies the stream used to print
556  /// banners (by default cout). If a null pointer is passed, banner
557  /// printout is suppressed. This affects all banners, including
558  /// those from plugins.
559  ///
560  /// Please note that if you distribute 3rd party code
561  /// that links with FastJet, that 3rd party code must not
562  /// use this call turn off the printing of FastJet banners
563  /// by default. This requirement reflects the spirit of
564  /// clause 2c of the GNU Public License (v2), under which
565  /// FastJet and its plugins are distributed.
566  static void set_fastjet_banner_stream(std::ostream * ostr) {_fastjet_banner_ostr = ostr;}
567  // [this line must be left as is to hide the doxygen comment]
568  /// \endcond
569 
570  /// returns a pointer to the stream to be used to print banners
571  /// (cout by default). This function is used by plugins to determine
572  /// where to direct their banners. Plugins should properly handle
573  /// the case where the pointer is null.
574  static std::ostream * fastjet_banner_stream() {return _fastjet_banner_ostr;}
575 
576 private:
577  /// \cond internal_doc
578 
579  /// contains the actual stream to use for banners
580  static std::ostream * _fastjet_banner_ostr;
581 
582  /// \endcond
583 
584 protected:
585 //DEP static JetAlgorithm _default_jet_algorithm;
586  JetDefinition _jet_def;
587 
588  /// transfer the vector<L> of input jets into our own vector<PseudoJet>
589  /// _jets (with some reserved space for future growth).
590  template<class L> void _transfer_input_jets(
591  const std::vector<L> & pseudojets);
592 
593  /// This is what is called to do all the initialisation and
594  /// then run the clustering (may be called by various constructors).
595  /// It assumes _jets contains the momenta to be clustered.
596  void _initialise_and_run (const JetDefinition & jet_def,
597  const bool & writeout_combinations);
598 
599  //// this performs the initialisation, minus the option-decanting
600  //// stage; for low multiplicity, initialising a few things in the
601  //// constructor, calling the decant_options_partial() and then this
602  //// is faster than going through _initialise_and_run.
603  void _initialise_and_run_no_decant();
604 
605 //DEP /// This is an alternative routine for initialising and running the
606 //DEP /// clustering, provided for legacy purposes. The jet finder is that
607 //DEP /// specified in the static member _default_jet_algorithm.
608 //DEP void _initialise_and_run (const double R,
609 //DEP const Strategy & strategy,
610 //DEP const bool & writeout_combinations);
611 
612  /// fills in the various member variables with "decanted" options from
613  /// the jet_definition and writeout_combinations variables
614  void _decant_options(const JetDefinition & jet_def,
615  const bool & writeout_combinations);
616 
617  /// assuming that the jet definition, writeout_combinations and
618  /// _structure_shared_ptr have been set (e.g. in an initialiser list
619  /// in the constructor), it handles the remaining decanting of
620  /// options.
621  void _decant_options_partial();
622 
623  /// fill out the history (and jet cross refs) related to the initial
624  /// set of jets (assumed already to have been "transferred"),
625  /// without any clustering
626  void _fill_initial_history();
627 
628  /// carry out the recombination between the jets numbered jet_i and
629  /// jet_j, at distance scale dij; return the index newjet_k of the
630  /// result of the recombination of i and j.
631  void _do_ij_recombination_step(const int jet_i, const int jet_j,
632  const double dij, int & newjet_k);
633 
634  /// carry out an recombination step in which _jets[jet_i] merges with
635  /// the beam,
636  void _do_iB_recombination_step(const int jet_i, const double diB);
637 
638  /// every time a jet is added internally during clustering, this
639  /// should be called to set the jet's structure shared ptr to point
640  /// to the CS (and the count of internally associated objects is
641  /// also updated). This should not be called outside construction of
642  /// a CS object.
643  void _set_structure_shared_ptr(PseudoJet & j);
644 
645  /// make sure that the CS's internal tally of the use count matches
646  /// that of the _structure_shared_ptr
647  void _update_structure_use_count();
648 
649  /// returns a suggestion for the best strategy to use on event
650  /// multiplicity, algorithm, R, etc.
651  Strategy _best_strategy() const;
652 
653  /// \if internal_doc
654  /// \class _Parabola
655  /// returns c*(a*R**2 + b*R + 1);
656  /// Written as a class in case we want to give names to different
657  /// parabolas
658  /// \endif
659  class _Parabola {
660  public:
661  _Parabola(double a, double b, double c) : _a(a), _b(b), _c(c) {}
662  inline double operator()(const double R) const {return _c*(_a*R*R + _b*R + 1);}
663  private:
664  double _a, _b, _c;
665  };
666 
667  /// \if internal_doc
668  /// \class _Line
669  /// operator()(R) returns a*R+b;
670  /// \endif
671  class _Line {
672  public:
673  _Line(double a, double b) : _a(a), _b(b) {}
674  inline double operator()(const double R) const {return _a*R + _b;}
675  private:
676  double _a, _b;
677  };
678 
679  /// This contains the physical PseudoJets; for each PseudoJet one
680  /// can find the corresponding position in the _history by looking
681  /// at _jets[i].cluster_hist_index().
682  std::vector<PseudoJet> _jets;
683 
684 
685  /// this vector will contain the branching history; for each stage,
686  /// _history[i].jetp_index indicates where to look in the _jets
687  /// vector to get the physical PseudoJet.
688  std::vector<history_element> _history;
689 
690  /// set subhist to be a set pointers to history entries corresponding to the
691  /// subjets of this jet; one stops going working down through the
692  /// subjets either when
693  /// - there is no further to go
694  /// - one has found maxjet entries
695  /// - max_dij_so_far <= dcut
696  /// By setting maxjet=0 one can use just dcut; by setting dcut<0
697  /// one can use jet maxjet
698  void get_subhist_set(std::set<const history_element*> & subhist,
699  const PseudoJet & jet, double dcut, int maxjet) const;
700 
701  bool _writeout_combinations;
702  int _initial_n;
703  double _Rparam, _R2, _invR2;
704  double _Qtot;
705  Strategy _strategy;
706  JetAlgorithm _jet_algorithm;
707 
708  SharedPtr<PseudoJetStructureBase> _structure_shared_ptr; //< will actually be of type ClusterSequenceStructure
709  int _structure_use_count_after_construction; //< info of use when CS handles its own memory
710  /// if true then the CS will delete itself when the last external
711  /// object referring to it disappears. It is mutable so as to ensure
712  /// that signal_imminent_self_deletion() [const] can make relevant
713  /// changes.
715 
716  private:
717 
718  bool _plugin_activated;
719  SharedPtr<Extras> _extras; // things the plugin might want to add
720 
721  void _really_dumb_cluster ();
722  void _delaunay_cluster ();
723  //void _simple_N2_cluster ();
724  template<class BJ> void _simple_N2_cluster ();
725  void _tiled_N2_cluster ();
726  void _faster_tiled_N2_cluster ();
727 
728  //
729  void _minheap_faster_tiled_N2_cluster();
730 
731  // things needed specifically for Cambridge with Chan's 2D closest
732  // pairs method
733  void _CP2DChan_cluster();
734  void _CP2DChan_cluster_2pi2R ();
735  void _CP2DChan_cluster_2piMultD ();
736  void _CP2DChan_limited_cluster(double D);
737  void _do_Cambridge_inclusive_jets();
738 
739  // NSqrtN method for C/A
740  void _fast_NsqrtN_cluster();
741 
742  void _add_step_to_history( //const int step_number,
743  const int parent1,
744  const int parent2, const int jetp_index,
745  const double dij);
746 
747  /// internal routine associated with the construction of the unique
748  /// history order (following children in the tree)
749  void _extract_tree_children(int pos, std::valarray<bool> &,
750  const std::valarray<int> &, std::vector<int> &) const;
751 
752  /// internal routine associated with the construction of the unique
753  /// history order (following parents in the tree)
754  void _extract_tree_parents (int pos, std::valarray<bool> &,
755  const std::valarray<int> &, std::vector<int> &) const;
756 
757 
758  // these will be useful shorthands in the Voronoi-based code
759  typedef std::pair<int,int> TwoVertices;
760  typedef std::pair<double,TwoVertices> DijEntry;
761  typedef std::multimap<double,TwoVertices> DistMap;
762 
763  /// currently used only in the Voronoi based code
764  void _add_ktdistance_to_map(const int ii,
765  DistMap & DijMap,
766  const DynamicNearestNeighbours * DNN);
767 
768 
769  /// will be set by default to be true for the first run
770  static bool _first_time;
771 
772  /// manage warnings related to exclusive jets access
773  static LimitedWarning _exclusive_warnings;
774 
775  /// the limited warning member for notification of user that
776  /// their requested strategy has been overridden (usually because
777  /// they have R>2pi and not all strategies work then)
778  static LimitedWarning _changed_strategy_warning;
779 
780  //----------------------------------------------------------------------
781  /// the fundamental structure which contains the minimal info about
782  /// a jet, as needed for our plain N^2 algorithm -- the idea is to
783  /// put all info that will be accessed N^2 times into an array of
784  /// BriefJets...
785  struct BriefJet {
786  double eta, phi, kt2, NN_dist;
787  BriefJet * NN;
788  int _jets_index;
789  };
790 
791  /// structure analogous to BriefJet, but with the extra information
792  /// needed for dealing with tiles
793  class TiledJet {
794  public:
795  double eta, phi, kt2, NN_dist;
796  TiledJet * NN, *previous, * next;
797  int _jets_index, tile_index, diJ_posn;
798  // routines that are useful in the minheap version of tiled
799  // clustering ("misuse" the otherwise unused diJ_posn, so as
800  // to indicate whether jets need to have their minheap entries
801  // updated).
802  inline void label_minheap_update_needed() {diJ_posn = 1;}
803  inline void label_minheap_update_done() {diJ_posn = 0;}
804  inline bool minheap_update_needed() const {return diJ_posn==1;}
805  };
806 
807  //-- some of the functions that follow are templates and will work
808  //as well for briefjet and tiled jets
809 
810  /// set the kinematic and labelling info for jeta so that it corresponds
811  /// to _jets[_jets_index]
812  template <class J> void _bj_set_jetinfo( J * const jet,
813  const int _jets_index) const;
814 
815  /// "remove" this jet, which implies updating links of neighbours and
816  /// perhaps modifying the tile structure
817  void _bj_remove_from_tiles( TiledJet * const jet) const;
818 
819  /// return the distance between two BriefJet objects
820  template <class J> double _bj_dist(const J * const jeta,
821  const J * const jetb) const;
822 
823  // return the diJ (multiplied by _R2) for this jet assuming its NN
824  // info is correct
825  template <class J> double _bj_diJ(const J * const jeta) const;
826 
827  /// for testing purposes only: if in the range head--tail-1 there is a
828  /// a jet which corresponds to hist_index in the history, then
829  /// return a pointer to that jet; otherwise return tail.
830  template <class J> inline J * _bj_of_hindex(
831  const int hist_index,
832  J * const head, J * const tail)
833  const {
834  J * res;
835  for(res = head; res<tail; res++) {
836  if (_jets[res->_jets_index].cluster_hist_index() == hist_index) {break;}
837  }
838  return res;
839  }
840 
841 
842  //-- remaining functions are different in various cases, so we
843  // will use templates but are not sure if they're useful...
844 
845  /// updates (only towards smaller distances) the NN for jeta without checking
846  /// whether in the process jeta itself might be a new NN of one of
847  /// the jets being scanned -- span the range head to tail-1 with
848  /// assumption that jeta is not contained in that range
849  template <class J> void _bj_set_NN_nocross(J * const jeta,
850  J * const head, const J * const tail) const;
851 
852  /// reset the NN for jeta and DO check whether in the process jeta
853  /// itself might be a new NN of one of the jets being scanned --
854  /// span the range head to tail-1 with assumption that jeta is not
855  /// contained in that range
856  template <class J> void _bj_set_NN_crosscheck(J * const jeta,
857  J * const head, const J * const tail) const;
858 
859 
860 
861  /// number of neighbours that a tile will have (rectangular geometry
862  /// gives 9 neighbours).
863  static const int n_tile_neighbours = 9;
864  //----------------------------------------------------------------------
865  /// The fundamental structures to be used for the tiled N^2 algorithm
866  /// (see CCN27-44 for some discussion of pattern of tiling)
867  struct Tile {
868  /// pointers to neighbouring tiles, including self
869  Tile * begin_tiles[n_tile_neighbours];
870  /// neighbouring tiles, excluding self
871  Tile ** surrounding_tiles;
872  /// half of neighbouring tiles, no self
873  Tile ** RH_tiles;
874  /// just beyond end of tiles
875  Tile ** end_tiles;
876  /// start of list of BriefJets contained in this tile
877  TiledJet * head;
878  /// sometimes useful to be able to tag a tile
879  bool tagged;
880  };
881  std::vector<Tile> _tiles;
882  double _tiles_eta_min, _tiles_eta_max;
883  double _tile_size_eta, _tile_size_phi;
884  int _n_tiles_phi,_tiles_ieta_min,_tiles_ieta_max;
885 
886  // reasonably robust return of tile index given ieta and iphi, in particular
887  // it works even if iphi is negative
888  inline int _tile_index (int ieta, int iphi) const {
889  // note that (-1)%n = -1 so that we have to add _n_tiles_phi
890  // before performing modulo operation
891  return (ieta-_tiles_ieta_min)*_n_tiles_phi
892  + (iphi+_n_tiles_phi) % _n_tiles_phi;
893  }
894 
895  // routines for tiled case, including some overloads of the plain
896  // BriefJet cases
897  int _tile_index(const double eta, const double phi) const;
898  void _tj_set_jetinfo ( TiledJet * const jet, const int _jets_index);
899  void _bj_remove_from_tiles(TiledJet * const jet);
900  void _initialise_tiles();
901  void _print_tiles(TiledJet * briefjets ) const;
902  void _add_neighbours_to_tile_union(const int tile_index,
903  std::vector<int> & tile_union, int & n_near_tiles) const;
904  void _add_untagged_neighbours_to_tile_union(const int tile_index,
905  std::vector<int> & tile_union, int & n_near_tiles);
906 
907  //----------------------------------------------------------------------
908  /// fundamental structure for e+e- clustering
909  struct EEBriefJet {
910  double NN_dist; // obligatorily present
911  double kt2; // obligatorily present == E^2 in general
912  EEBriefJet * NN; // must be present too
913  int _jets_index; // must also be present!
914  //...........................................................
915  double nx, ny, nz; // our internal storage for fast distance calcs
916  };
917 
918  /// to help instantiation (fj 2.4.0; did not quite work on gcc 33 and os x 10.3?)
919  //void _dummy_N2_cluster_instantiation();
920 
921 
922  /// to avoid issues with template instantiation (OS X 10.3, gcc 3.3)
923  void _simple_N2_cluster_BriefJet();
924  /// to avoid issues with template instantiation (OS X 10.3, gcc 3.3)
925  void _simple_N2_cluster_EEBriefJet();
926 };
927 
928 
929 //**********************************************************************
930 //************** START OF INLINE MATERIAL ******************
931 //**********************************************************************
932 
933 
934 //----------------------------------------------------------------------
935 // Transfer the initial jets into our internal structure
936 template<class L> void ClusterSequence::_transfer_input_jets(
937  const std::vector<L> & pseudojets) {
938 
939  // this will ensure that we can point to jets without difficulties
940  // arising.
941  _jets.reserve(pseudojets.size()*2);
942 
943  // insert initial jets this way so that any type L that can be
944  // converted to a pseudojet will work fine (basically PseudoJet
945  // and any type that has [] subscript access to the momentum
946  // components, such as CLHEP HepLorentzVector).
947  for (unsigned int i = 0; i < pseudojets.size(); i++) {
948  _jets.push_back(pseudojets[i]);}
949 
950 }
951 
952 // //----------------------------------------------------------------------
953 // // initialise from some generic type... Has to be made available
954 // // here in order for it the template aspect of it to work...
955 // template<class L> ClusterSequence::ClusterSequence (
956 // const std::vector<L> & pseudojets,
957 // const double R,
958 // const Strategy & strategy,
959 // const bool & writeout_combinations) {
960 //
961 // // transfer the initial jets (type L) into our own array
962 // _transfer_input_jets(pseudojets);
963 //
964 // // run the clustering
965 // _initialise_and_run(R,strategy,writeout_combinations);
966 // }
967 
968 
969 //----------------------------------------------------------------------
970 /// constructor of a jet-clustering sequence from a vector of
971 /// four-momenta, with the jet definition specified by jet_def
972 template<class L> ClusterSequence::ClusterSequence (
973  const std::vector<L> & pseudojets,
974  const JetDefinition & jet_def_in,
975  const bool & writeout_combinations) :
976  _jet_def(jet_def_in), _writeout_combinations(writeout_combinations),
977  _structure_shared_ptr(new ClusterSequenceStructure(this))
978 {
979 
980  // transfer the initial jets (type L) into our own array
981  _transfer_input_jets(pseudojets);
982 
983  // transfer the remaining options
985 
986  // run the clustering
987  _initialise_and_run_no_decant();
988 }
989 
990 
991 inline const std::vector<PseudoJet> & ClusterSequence::jets () const {
992  return _jets;
993 }
994 
995 inline const std::vector<ClusterSequence::history_element> & ClusterSequence::history () const {
996  return _history;
997 }
998 
999 inline unsigned int ClusterSequence::n_particles() const {return _initial_n;}
1000 
1001 //----------------------------------------------------------------------
1002 // implementation of JetDefinition::operator() is here to avoid nasty
1003 // issues of order of implementations and includes
1004 #ifndef __CINT__
1005 template<class L>
1006 std::vector<PseudoJet> JetDefinition::operator()(const std::vector<L> & particles) const {
1007  // create a new cluster sequence
1008  ClusterSequence * cs = new ClusterSequence(particles, *this);
1009 
1010  // get the jets, and sort them according to whether the algorithm
1011  // is spherical or not
1012  std::vector<PseudoJet> jets;
1013  if (is_spherical()) {
1014  jets = sorted_by_E(cs->inclusive_jets());
1015  } else {
1016  jets = sorted_by_pt(cs->inclusive_jets());
1017  }
1018 
1019  // make sure the ClusterSequence gets deleted once it's no longer
1020  // needed
1021  if (jets.size() != 0) {
1023  } else {
1024  delete cs;
1025  }
1026 
1027  return jets;
1028 }
1029 #endif // __CINT__
1030 
1031 
1032 //----------------------------------------------------------------------
1033 template <class J> inline void ClusterSequence::_bj_set_jetinfo(
1034  J * const jetA, const int _jets_index) const {
1035  jetA->eta = _jets[_jets_index].rap();
1036  jetA->phi = _jets[_jets_index].phi_02pi();
1037  jetA->kt2 = jet_scale_for_algorithm(_jets[_jets_index]);
1038  jetA->_jets_index = _jets_index;
1039  // initialise NN info as well
1040  jetA->NN_dist = _R2;
1041  jetA->NN = NULL;
1042 }
1043 
1044 
1045 
1046 
1047 //----------------------------------------------------------------------
1048 template <class J> inline double ClusterSequence::_bj_dist(
1049  const J * const jetA, const J * const jetB) const {
1050  //#define FASTJET_NEW_DELTA_PHI
1051 #ifndef FASTJET_NEW_DELTA_PHI
1052  //GPS+MC old version of Delta phi calculation
1053  double dphi = std::abs(jetA->phi - jetB->phi);
1054  double deta = (jetA->eta - jetB->eta);
1055  if (dphi > pi) {dphi = twopi - dphi;}
1056 #else
1057  //GPS+MC testing for 2015-02-faster-deltaR2
1058  double dphi = pi-std::abs(pi-std::abs(jetA->phi - jetB->phi));
1059  double deta = (jetA->eta - jetB->eta);
1060 #endif
1061  return dphi*dphi + deta*deta;
1062 }
1063 
1064 //----------------------------------------------------------------------
1065 template <class J> inline double ClusterSequence::_bj_diJ(const J * const jet) const {
1066  double kt2 = jet->kt2;
1067  if (jet->NN != NULL) {if (jet->NN->kt2 < kt2) {kt2 = jet->NN->kt2;}}
1068  return jet->NN_dist * kt2;
1069 }
1070 
1071 
1072 //----------------------------------------------------------------------
1073 // set the NN for jet without checking whether in the process you might
1074 // have discovered a new nearest neighbour for another jet
1075 template <class J> inline void ClusterSequence::_bj_set_NN_nocross(
1076  J * const jet, J * const head, const J * const tail) const {
1077  double NN_dist = _R2;
1078  J * NN = NULL;
1079  if (head < jet) {
1080  for (J * jetB = head; jetB != jet; jetB++) {
1081  double dist = _bj_dist(jet,jetB);
1082  if (dist < NN_dist) {
1083  NN_dist = dist;
1084  NN = jetB;
1085  }
1086  }
1087  }
1088  if (tail > jet) {
1089  for (J * jetB = jet+1; jetB != tail; jetB++) {
1090  double dist = _bj_dist(jet,jetB);
1091  if (dist < NN_dist) {
1092  NN_dist = dist;
1093  NN = jetB;
1094  }
1095  }
1096  }
1097  jet->NN = NN;
1098  jet->NN_dist = NN_dist;
1099 }
1100 
1101 
1102 //----------------------------------------------------------------------
1103 template <class J> inline void ClusterSequence::_bj_set_NN_crosscheck(J * const jet,
1104  J * const head, const J * const tail) const {
1105  double NN_dist = _R2;
1106  J * NN = NULL;
1107  for (J * jetB = head; jetB != tail; jetB++) {
1108  double dist = _bj_dist(jet,jetB);
1109  if (dist < NN_dist) {
1110  NN_dist = dist;
1111  NN = jetB;
1112  }
1113  if (dist < jetB->NN_dist) {
1114  jetB->NN_dist = dist;
1115  jetB->NN = jet;
1116  }
1117  }
1118  jet->NN = NN;
1119  jet->NN_dist = NN_dist;
1120 }
1121 
1122 FASTJET_END_NAMESPACE
1123 
1124 #endif // __FASTJET_CLUSTERSEQUENCE_HH__
const JetDefinition & jet_def() const
return a reference to the jet definition
ClusterSequenceStructure StructureType
the structure type associated with a jet belonging to 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...
void plugin_associate_extras(Extras *extras_in)
the plugin can associate some extra information with the ClusterSequence object by calling this funct...
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...
vector< PseudoJet > sorted_by_pt(const vector< PseudoJet > &jets)
return a vector of jets sorted into decreasing kt2
Definition: PseudoJet.cc:774
const Extras * extras() const
returns a pointer to the extras object (may be null)
deals with clustering
const std::vector< history_element > & history() const
allow the user to access the raw internal history.
void plugin_simple_N2_cluster()
allows a plugin to run a templated clustering (nearest-neighbour heuristic)
ClusterSequence()
default constructor
ClusterSequence(const ClusterSequence &cs)
copy constructor for a ClusterSequence
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.
int jetp_index
index in _history where the current jet is recombined with another jet to form its child...
base class to store extra information that plugins may provide
double Q() const
returns the sum of all energies in the event (relevant mainly for e+e-)
double Q2() const
return Q()^2
Contains any information related to the clustering that should be directly accessible to PseudoJet...
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], with the specified dij, and return the index (newjet_k) allocated to the new jet, whose momentum is assumed to be the 4-vector sum of that of jet_i and jet_j
std::vector< PseudoJet > _jets
This contains the physical PseudoJets; for each PseudoJet one can find the corresponding position in ...
std::vector< history_element > _history
this vector will contain the branching history; for each stage, _history[i].jetp_index indicates wher...
bool plugin_activated() const
returns true when the plugin is allowed to run the show.
class to provide facilities for giving warnings up to some maximum number of times and to provide glo...
int parent2
index in _history where first parent of this jet was created (InexistentParent if this jet is an orig...
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...
double dij
index in the _jets vector where we will find the
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...
std::vector< PseudoJet > exclusive_jets_ycut(double ycut) const
the exclusive jets obtained at the given ycut
double exclusive_ymerge_max(int njets) const
same as exclusive_dmerge_max, but normalised to squared total energy
vector< PseudoJet > sorted_by_E(const vector< PseudoJet > &jets)
return a vector of jets sorted into decreasing energy
Definition: PseudoJet.cc:790
Strategy strategy_used() const
return the enum value of the strategy used to cluster the event
int child
index in _history where second parent of this jet was created (InexistentParent if this jet is an ori...
const SharedPtr< PseudoJetStructureBase > & structure_shared_ptr() const
retrieve a shared pointer to the wrapper to this ClusterSequence
an implementation of C++0x shared pointers (or boost&#39;s)
Definition: SharedPtr.hh:121
std::string strategy_string() const
return the name of the strategy used to cluster the event
a single element in the clustering history
int n_exclusive_jets_ycut(double ycut) const
the number of exclusive jets at the given ycut
bool will_delete_self_when_unused() const
return true if the object has been told to delete itself when unused
static std::ostream * fastjet_banner_stream()
returns a pointer to the stream to be used to print banners (cout by default).
void _decant_options_partial()
assuming that the jet definition, writeout_combinations and _structure_shared_ptr have been set (e...
JetAlgorithm
the various families of jet-clustering algorithm
void delete_self_when_unused()
by calling this routine you tell the ClusterSequence to delete itself when all the Pseudojets associa...
Strategy
the various options for the algorithmic strategy to adopt in clustering events with kt and cambridge ...
bool _deletes_self_when_unused
if true then the CS will delete itself when the last external object referring to it disappears...
Class to contain pseudojets, including minimal information of use to jet-clustering routines...
Definition: PseudoJet.hh:67
double max_dij_so_far
the distance corresponding to the recombination at this stage of the clustering.
const std::vector< PseudoJet > & jets() const
allow the user to access the internally stored _jets() array, which contains both the initial particl...
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...
class that is intended to hold a full definition of the jet clusterer
unsigned int n_particles() const
returns the number of particles that were provided to the clustering algorithm (helps the user find t...