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