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