1 | //STARTHEADER |
---|
2 | // $Id: ClusterSequence.hh 2806 2011-12-01 17:21: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/internal/DynamicNearestNeighbours.hh" |
---|
36 | #include "fastjet/PseudoJet.hh" |
---|
37 | #include<memory> |
---|
38 | #include<cassert> |
---|
39 | #include<iostream> |
---|
40 | #include<string> |
---|
41 | #include<set> |
---|
42 | #include<cmath> // needed to get double std::abs(double) |
---|
43 | #include "fastjet/Error.hh" |
---|
44 | #include "fastjet/JetDefinition.hh" |
---|
45 | #include "fastjet/SharedPtr.hh" |
---|
46 | #include "fastjet/LimitedWarning.hh" |
---|
47 | #include "fastjet/FunctionOfPseudoJet.hh" |
---|
48 | #include "fastjet/ClusterSequenceStructure.hh" |
---|
49 | |
---|
50 | FASTJET_BEGIN_NAMESPACE // defined in fastjet/internal/base.hh |
---|
51 | |
---|
52 | |
---|
53 | // forward declaration |
---|
54 | class ClusterSequenceStructure; |
---|
55 | |
---|
56 | /// @ingroup basic_classes |
---|
57 | /// \class ClusterSequence |
---|
58 | /// deals with clustering |
---|
59 | class ClusterSequence { |
---|
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). |
---|
404 | struct history_element{ |
---|
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 |
---|
530 | const SharedPtr<PseudoJetStructureBase> & structure_shared_ptr() const{ |
---|
531 | return _structure_shared_ptr; |
---|
532 | } |
---|
533 | |
---|
534 | /// the structure type associated with a jet belonging to a ClusterSequence |
---|
535 | typedef ClusterSequenceStructure StructureType; |
---|
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. |
---|
678 | mutable bool _deletes_self_when_unused; |
---|
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 |
---|
952 | _decant_options_partial(); |
---|
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__ |
---|