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