Fork me on GitHub

source: git/external/fastjet/ClusterSequence.hh@ 0e6b14e

ImprovedOutputFile Timing
Last change on this file since 0e6b14e was b7b836a, checked in by Pavel Demin <pavel-demin@…>, 7 years ago

update FastJet library to 3.3.1 and FastJet Contrib library to 1.036

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