Fork me on GitHub

source: svn/trunk/Utilities/Fastjet/include/fastjet/ClusterSequence.hh@ 480

Last change on this file since 480 was 11, checked in by severine ovyn, 16 years ago

Fastjet added; CDFCones directory has been changed

File size: 28.8 KB
RevLine 
[11]1//STARTHEADER
2// $Id: ClusterSequence.hh,v 1.1 2008-11-06 14:32:07 ovyn Exp $
3//
4// Copyright (c) 2005-2006, Matteo Cacciari and Gavin Salam
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, write to the Free Software
26// Foundation, Inc.:
27// 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
28//----------------------------------------------------------------------
29//ENDHEADER
30
31
32//----------------------------------------------------------------------
33// here's where we put the main page for fastjet (as explained in the
34// Doxygen faq)
35//......................................................................
36/*! \mainpage FastJet code documentation
37 *
38 * These pages provide automatically generated documentation for the
39 * FastJet package.
40 *
41 * For further information and normal documentation, see the main <a
42 * href="http://www.lpthe.jussieu.fr/~salam/fastjet">FastJet</a> page.
43 */
44//----------------------------------------------------------------------
45
46#ifndef __FASTJET_CLUSTERSEQUENCE_HH__
47#define __FASTJET_CLUSTERSEQUENCE_HH__
48
49#include<vector>
50#include<map>
51#include "Utilities/Fastjet/include/fastjet/internal/DynamicNearestNeighbours.hh"
52#include "Utilities/Fastjet/include/fastjet/PseudoJet.hh"
53#include<memory>
54#include<cassert>
55#include<iostream>
56#include<string>
57#include<cmath> // needed to get double std::abs(double)
58#include "Utilities/Fastjet/include/fastjet/Error.hh"
59#include "Utilities/Fastjet/include/fastjet/JetDefinition.hh"
60
61FASTJET_BEGIN_NAMESPACE // defined in fastjet/internal/base.hh
62
63
64/// deals with clustering
65class ClusterSequence {
66
67
68 public:
69
70 /// default constructor
71 ClusterSequence () {}
72
73 /// create a clustersequence starting from the supplied set
74 /// of pseudojets and clustering them with the long-invariant
75 /// kt algorithm (E-scheme recombination) with the supplied
76 /// value for R.
77 ///
78 /// If strategy=DumbN3 a very stupid N^3 algorithm is used for the
79 /// clustering; otherwise strategy = NlnN* uses cylinders algorithms
80 /// with some number of pi coverage. If writeout_combinations=true a
81 /// summary of the recombination sequence is written out
82 template<class L> ClusterSequence (const std::vector<L> & pseudojets,
83 const double & R = 1.0,
84 const Strategy & strategy = Best,
85 const bool & writeout_combinations = false);
86
87
88 /// create a clustersequence starting from the supplied set
89 /// of pseudojets and clustering them with jet definition specified
90 /// by jet_def (which also specifies the clustering strategy)
91 template<class L> ClusterSequence (
92 const std::vector<L> & pseudojets,
93 const JetDefinition & jet_def,
94 const bool & writeout_combinations = false);
95
96
97 // NB: in the routines that follow, for extracting lists of jets, a
98 // list structure might be more efficient, if sometimes a little
99 // more awkward to use (at least for old fortran hands).
100
101 /// return a vector of all jets (in the sense of the inclusive
102 /// algorithm) with pt >= ptmin. Time taken should be of the order
103 /// of the number of jets returned.
104 std::vector<PseudoJet> inclusive_jets (const double & ptmin = 0.0) const;
105
106 /// return the number of jets (in the sense of the exclusive
107 /// algorithm) that would be obtained when running the algorithm
108 /// with the given dcut.
109 int n_exclusive_jets (const double & dcut) const;
110
111 /// return a vector of all jets (in the sense of the exclusive
112 /// algorithm) that would be obtained when running the algorithm
113 /// with the given dcut.
114 std::vector<PseudoJet> exclusive_jets (const double & dcut) const;
115
116 /// return a vector of all jets when the event is clustered (in the
117 /// exclusive sense) to exactly njets.
118 std::vector<PseudoJet> exclusive_jets (const int & njets) const;
119
120 /// return the dmin corresponding to the recombination that went from
121 /// n+1 to n jets (sometimes known as d_{n n+1}).
122 double exclusive_dmerge (const int & njets) const;
123
124 /// return the maximum of the dmin encountered during all recombinations
125 /// up to the one that led to an n-jet final state; identical to
126 /// exclusive_dmerge, except in cases where the dmin do not increase
127 /// monotonically.
128 double exclusive_dmerge_max (const int & njets) const;
129
130 /// returns true iff the object is included in the jet.
131 ///
132 /// NB: this is only sensible if the object is already registered
133 /// within the cluster sequence, so you cannot use it with an input
134 /// particle to the CS (since the particle won't have the history
135 /// index set properly).
136 ///
137 /// For nice clustering structures it should run in O(ln(N)) time
138 /// but in worst cases (certain cone plugins) it can take O(n) time,
139 /// where n is the number of particles in the jet.
140 bool object_in_jet(const PseudoJet & object, const PseudoJet & jet) const;
141
142 /// if the jet has parents in the clustering, it returns true
143 /// and sets parent1 and parent2 equal to them.
144 ///
145 /// if it has no parents it returns false and sets parent1 and
146 /// parent2 to zero
147 bool has_parents(const PseudoJet & jet, PseudoJet & parent1,
148 PseudoJet & parent2) const;
149
150 /// if the jet has a child then return true and give the child jet
151 /// otherwise return false and set the child to zero
152 bool has_child(const PseudoJet & jet, PseudoJet & child) const;
153
154 /// Version of has_child that sets a pointer to the child if the child
155 /// exists;
156 bool has_child(const PseudoJet & jet, const PseudoJet * & childp) const;
157
158 /// if this jet has a child (and so a partner) return true
159 /// and give the partner, otherwise return false and set the
160 /// partner to zero
161 bool has_partner(const PseudoJet & jet, PseudoJet & partner) const;
162
163 /// return a vector of the particles that make up jet
164 std::vector<PseudoJet> constituents (const PseudoJet & jet) const;
165
166
167 /// output the supplied vector of jets in a format that can be read
168 /// by an appropriate root script; the format is:
169 /// jet-n jet-px jet-py jet-pz jet-E
170 /// particle-n particle-rap particle-phi particle-pt
171 /// particle-n particle-rap particle-phi particle-pt
172 /// ...
173 /// #END
174 /// ... [i.e. above repeated]
175 void print_jets_for_root(const std::vector<PseudoJet> & jets,
176 std::ostream & ostr = std::cout) const;
177
178// Not yet. Perhaps in a future release.
179// /// print out all inclusive jets with pt > ptmin
180// virtual void print_jets (const double & ptmin=0.0) const;
181
182 /// add on to subjet_vector the subjets of jet (for internal use mainly)
183 void add_constituents (const PseudoJet & jet,
184 std::vector<PseudoJet> & subjet_vector) const;
185
186 /// return the enum value of the strategy used to cluster the event
187 inline Strategy strategy_used () const {return _strategy;}
188 std::string strategy_string () const;
189
190 /// return a reference to the jet definition
191 const JetDefinition & jet_def() const {return _jet_def;}
192
193 /// returns the scale associated with a jet as required for this
194 /// clustering algorithm (kt^2 for the kt-algorithm, 1 for the
195 /// Cambridge algorithm). [May become virtual at some point]
196 double jet_scale_for_algorithm(const PseudoJet & jet) const;
197
198 //----- next follow functions designed specifically for plugins, which
199 // may only be called when plugin_activated() returns true
200
201 /// record the fact that there has been a recombination between
202 /// jets()[jet_i] and jets()[jet_k], with the specified dij, and
203 /// return the index (newjet_k) allocated to the new jet, whose
204 /// momentum is assumed to be the 4-vector sum of that of jet_i and
205 /// jet_j
206 void plugin_record_ij_recombination(int jet_i, int jet_j, double dij,
207 int & newjet_k) {
208 assert(plugin_activated());
209 _do_ij_recombination_step(jet_i, jet_j, dij, newjet_k);
210 }
211
212 /// as for the simpler variant of plugin_record_ij_recombination,
213 /// except that the new jet is attributed the momentum and
214 /// user_index of newjet
215 void plugin_record_ij_recombination(int jet_i, int jet_j, double dij,
216 const PseudoJet & newjet,
217 int & newjet_k);
218
219 /// record the fact that there has been a recombination between
220 /// jets()[jet_i] and the beam, with the specified diB; when looking
221 /// for inclusive jets, any iB recombination will returned to the user
222 /// as a jet.
223 void plugin_record_iB_recombination(int jet_i, double diB) {
224 assert(plugin_activated());
225 _do_iB_recombination_step(jet_i, diB);
226 }
227
228 /// a class intended to serve as a base in case a plugin needs to
229 /// associate extra information with a ClusterSequence (see
230 /// SISConePlugin.* for an example).
231 class Extras {
232 public:
233 virtual ~Extras() {}
234 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";}
235 };
236
237 /// the plugin can associated some extra information with the
238 /// ClusterSequence object by calling this function
239 inline void plugin_associate_extras(std::auto_ptr<Extras> extras_in) {
240 _extras = extras_in;
241 }
242
243 /// returns true when the plugin is allowed to run the show.
244 inline bool plugin_activated() const {return _plugin_activated;}
245
246 /// returns a pointer to the extras object (may be null)
247 const Extras * extras() const {return _extras.get();}
248
249public:
250 /// set the default (static) jet finder across all current and future
251 /// ClusterSequence objects -- deprecated and obsolescent (i.e. may be
252 /// suppressed in a future release).
253 static void set_jet_algorithm (JetAlgorithm jet_algorithm) {_default_jet_algorithm = jet_algorithm;}
254 /// same as above for backward compatibility
255 static void set_jet_finder (JetAlgorithm jet_algorithm) {_default_jet_algorithm = jet_algorithm;}
256
257
258 /// a single element in the clustering history (see vector _history
259 /// below).
260 struct history_element{
261 int parent1; /// index in _history where first parent of this jet
262 /// was created (InexistentParent if this jet is an
263 /// original particle)
264
265 int parent2; /// index in _history where second parent of this jet
266 /// was created (InexistentParent if this jet is an
267 /// original particle); BeamJet if this history entry
268 /// just labels the fact that the jet has recombined
269 /// with the beam)
270
271 int child; /// index in _history where the current jet is
272 /// recombined with another jet to form its child. It
273 /// is Invalid if this jet does not further
274 /// recombine.
275
276 int jetp_index; /// index in the _jets vector where we will find the
277 /// PseudoJet object corresponding to this jet
278 /// (i.e. the jet created at this entry of the
279 /// history). NB: if this element of the history
280 /// corresponds to a beam recombination, then
281 /// jetp_index=Invalid.
282
283 double dij; /// the distance corresponding to the recombination
284 /// at this stage of the clustering.
285
286 double max_dij_so_far; /// the largest recombination distance seen
287 /// so far in the clustering history.
288 };
289
290 enum JetType {Invalid=-3, InexistentParent = -2, BeamJet = -1};
291
292 /// allow the user to access the jets in this raw manner (needed
293 /// because we don't seem to be able to access protected elements of
294 /// the class for an object that is not "this" (at least in case where
295 /// "this" is of a slightly different kind from the object, both
296 /// derived from ClusterSequence).
297 const std::vector<PseudoJet> & jets() const;
298
299 /// allow the user to access the history in this raw manner (see
300 /// above for motivation).
301 const std::vector<history_element> & history() const;
302
303 /// returns the number of particles that were provided to the
304 /// clustering algorithm (helps the user find their way around the
305 /// history and jets objects if they weren't paying attention
306 /// beforehand).
307 unsigned int n_particles() const;
308
309 /// returns a vector of size n_particles() which indicates, for
310 /// each of the initial particles (in the order in which they were
311 /// supplied), which of the supplied jets it belongs to; if it does
312 /// not belong to any of the supplied jets, the index is set to -1;
313 std::vector<int> particle_jet_indices(const std::vector<PseudoJet> &) const;
314
315 /// routine that returns an order in which to read the history
316 /// such that clusterings that lead to identical jet compositions
317 /// but different histories (because of degeneracies in the
318 /// clustering order) will have matching constituents for each
319 /// matching entry in the unique_history_order.
320 ///
321 /// The order has the property that an entry's parents will always
322 /// appear prior to that entry itself.
323 ///
324 /// Roughly speaking the order is such that we first provide all
325 /// steps that lead to the final jet containing particle 1; then we
326 /// have the steps that lead to reconstruction of the jet containing
327 /// the next-lowest-numbered unclustered particle, etc...
328 /// [see GPS CCN28-12 for more info -- of course a full explanation
329 /// here would be better...]
330 std::vector<int> unique_history_order() const;
331
332 /// return the set of particles that have not been clustered. For
333 /// kt and cam/aachen algorithms this should always be null, but for
334 /// cone type algorithms it can be non-null;
335 std::vector<PseudoJet> unclustered_particles() const;
336
337 /// transfer the sequence contained in other_seq into our own;
338 /// any plugin "extras" contained in the from_seq will be lost
339 /// from there.
340 void transfer_from_sequence(ClusterSequence & from_seq);
341
342protected:
343 static JetAlgorithm _default_jet_algorithm;
344 JetDefinition _jet_def;
345
346 /// returns true if the jet has a history index contained within
347 /// the range of this CS
348 bool _potentially_valid(const PseudoJet & jet) const {
349 return jet.cluster_hist_index() >= 0
350 && jet.cluster_hist_index() < int(_history.size());
351 }
352
353 /// transfer the vector<L> of input jets into our own vector<PseudoJet>
354 /// _jets (with some reserved space for future growth).
355 template<class L> void _transfer_input_jets(
356 const std::vector<L> & pseudojets);
357
358 /// This is the routine that will do all the initialisation and
359 /// then run the clustering (may be called by various constructors).
360 /// It assumes _jets contains the momenta to be clustered.
361 void _initialise_and_run (const JetDefinition & jet_def,
362 const bool & writeout_combinations);
363
364 /// This is an alternative routine for initialising and running the
365 /// clustering, provided for legacy purposes. The jet finder is that
366 /// specified in the static member _default_jet_algorithm.
367 void _initialise_and_run (const double & R,
368 const Strategy & strategy,
369 const bool & writeout_combinations);
370
371 /// fills in the various member variables with "decanted" options from
372 /// the jet_definition and writeout_combinations variables
373 void _decant_options(const JetDefinition & jet_def,
374 const bool & writeout_combinations);
375
376 /// fill out the history (and jet cross refs) related to the initial
377 /// set of jets (assumed already to have been "transferred"),
378 /// without any clustering
379 void _fill_initial_history();
380
381 /// carry out the recombination between the jets numbered jet_i and
382 /// jet_j, at distance scale dij; return the index newjet_k of the
383 /// result of the recombination of i and j.
384 void _do_ij_recombination_step(const int & jet_i, const int & jet_j,
385 const double & dij, int & newjet_k);
386
387 /// carry out an recombination step in which _jets[jet_i] merges with
388 /// the beam,
389 void _do_iB_recombination_step(const int & jet_i, const double & diB);
390
391
392 /// This contains the physical PseudoJets; for each PseudoJet one
393 /// can find the corresponding position in the _history by looking
394 /// at _jets[i].cluster_hist_index().
395 std::vector<PseudoJet> _jets;
396
397
398 /// this vector will contain the branching history; for each stage,
399 /// _history[i].jetp_index indicates where to look in the _jets
400 /// vector to get the physical PseudoJet.
401 std::vector<history_element> _history;
402
403 bool _writeout_combinations;
404 int _initial_n;
405 double _Rparam, _R2, _invR2;
406 Strategy _strategy;
407 JetAlgorithm _jet_algorithm;
408
409 private:
410
411 bool _plugin_activated;
412 std::auto_ptr<Extras> _extras; // things the plugin might want to add
413
414 void _really_dumb_cluster ();
415 void _delaunay_cluster ();
416 void _simple_N2_cluster ();
417 void _tiled_N2_cluster ();
418 void _faster_tiled_N2_cluster ();
419
420 //
421 void _minheap_faster_tiled_N2_cluster();
422
423 // things needed specifically for Cambridge with Chan's 2D closest
424 // pairs method
425 void _CP2DChan_cluster();
426 void _CP2DChan_cluster_2pi2R ();
427 void _CP2DChan_cluster_2piMultD ();
428 void _CP2DChan_limited_cluster(double D);
429 void _do_Cambridge_inclusive_jets();
430
431 void _add_step_to_history(const int & step_number, const int & parent1,
432 const int & parent2, const int & jetp_index,
433 const double & dij);
434
435 /// internal routine associated with the construction of the unique
436 /// history order (following children in the tree)
437 void _extract_tree_children(int pos, std::valarray<bool> &,
438 const std::valarray<int> &, std::vector<int> &) const;
439
440 /// internal routine associated with the construction of the unique
441 /// history order (following parents in the tree)
442 void _extract_tree_parents (int pos, std::valarray<bool> &,
443 const std::valarray<int> &, std::vector<int> &) const;
444
445
446 // these will be useful shorthands in the Voronoi-based code
447 typedef std::pair<int,int> TwoVertices;
448 typedef std::pair<double,TwoVertices> DijEntry;
449 typedef std::multimap<double,TwoVertices> DistMap;
450
451 /// currently used only in the Voronoi based code
452 void _add_ktdistance_to_map(const int & ii,
453 DistMap & DijMap,
454 const DynamicNearestNeighbours * DNN);
455
456 /// for making sure the user knows what it is they're running...
457 void _print_banner();
458 /// will be set by default to be true for the first run
459 static bool _first_time;
460
461 /// record the number of warnings provided about the exclusive
462 /// algorithm -- so that we don't print it out more than a few
463 /// times.
464 static int _n_exclusive_warnings;
465
466 //----------------------------------------------------------------------
467 /// the fundamental structure which contains the minimal info about
468 /// a jet, as needed for our plain N^2 algorithm -- the idea is to
469 /// put all info that will be accessed N^2 times into an array of
470 /// BriefJets...
471 struct BriefJet {
472 double eta, phi, kt2, NN_dist;
473 BriefJet * NN;
474 int _jets_index;
475 };
476 /// structure analogous to BriefJet, but with the extra information
477 /// needed for dealing with tiles
478 class TiledJet {
479 public:
480 double eta, phi, kt2, NN_dist;
481 TiledJet * NN, *previous, * next;
482 int _jets_index, tile_index, diJ_posn;
483 // routines that are useful in the minheap version of tiled
484 // clustering ("misuse" the otherwise unused diJ_posn, so as
485 // to indicate whether jets need to have their minheap entries
486 // updated).
487 inline void label_minheap_update_needed() {diJ_posn = 1;}
488 inline void label_minheap_update_done() {diJ_posn = 0;}
489 inline bool minheap_update_needed() const {return diJ_posn==1;}
490 };
491
492 //-- some of the functions that follow are templates and will work
493 //as well for briefjet and tiled jets
494
495 /// set the kinematic and labelling info for jeta so that it corresponds
496 /// to _jets[_jets_index]
497 template <class J> void _bj_set_jetinfo( J * const jet,
498 const int _jets_index) const;
499
500 /// "remove" this jet, which implies updating links of neighbours and
501 /// perhaps modifying the tile structure
502 void _bj_remove_from_tiles( TiledJet * const jet) const;
503
504 /// return the distance between two BriefJet objects
505 template <class J> double _bj_dist(const J * const jeta,
506 const J * const jetb) const;
507
508 // return the diJ (multiplied by _R2) for this jet assuming its NN
509 // info is correct
510 template <class J> double _bj_diJ(const J * const jeta) const;
511
512 /// for testing purposes only: if in the range head--tail-1 there is a
513 /// a jet which corresponds to hist_index in the history, then
514 /// return a pointer to that jet; otherwise return tail.
515 template <class J> inline J * _bj_of_hindex(
516 const int hist_index,
517 J * const head, J * const tail)
518 const {
519 J * res;
520 for(res = head; res<tail; res++) {
521 if (_jets[res->_jets_index].cluster_hist_index() == hist_index) {break;}
522 }
523 return res;
524 }
525
526
527 //-- remaining functions are different in various cases, so we
528 // will use templates but are not sure if they're useful...
529
530 /// updates (only towards smaller distances) the NN for jeta without checking
531 /// whether in the process jeta itself might be a new NN of one of
532 /// the jets being scanned -- span the range head to tail-1 with
533 /// assumption that jeta is not contained in that range
534 template <class J> void _bj_set_NN_nocross(J * const jeta,
535 J * const head, const J * const tail) const;
536
537 /// reset the NN for jeta and DO check whether in the process jeta
538 /// itself might be a new NN of one of the jets being scanned --
539 /// span the range head to tail-1 with assumption that jeta is not
540 /// contained in that range
541 template <class J> void _bj_set_NN_crosscheck(J * const jeta,
542 J * const head, const J * const tail) const;
543
544
545
546 /// number of neighbours that a tile will have (rectangular geometry
547 /// gives 9 neighbours).
548 static const int n_tile_neighbours = 9;
549 //----------------------------------------------------------------------
550 /// The fundamental structures to be used for the tiled N^2 algorithm
551 /// (see CCN27-44 for some discussion of pattern of tiling)
552 struct Tile {
553 /// pointers to neighbouring tiles, including self
554 Tile * begin_tiles[n_tile_neighbours];
555 /// neighbouring tiles, excluding self
556 Tile ** surrounding_tiles;
557 /// half of neighbouring tiles, no self
558 Tile ** RH_tiles;
559 /// just beyond end of tiles
560 Tile ** end_tiles;
561 /// start of list of BriefJets contained in this tile
562 TiledJet * head;
563 /// sometimes useful to be able to tag a tile
564 bool tagged;
565 };
566 std::vector<Tile> _tiles;
567 double _tiles_eta_min, _tiles_eta_max;
568 double _tile_size_eta, _tile_size_phi;
569 int _n_tiles_phi,_tiles_ieta_min,_tiles_ieta_max;
570
571 // reasonably robust return of tile index given ieta and iphi, in particular
572 // it works even if iphi is negative
573 inline int _tile_index (int ieta, int iphi) const {
574 // note that (-1)%n = -1 so that we have to add _n_tiles_phi
575 // before performing modulo operation
576 return (ieta-_tiles_ieta_min)*_n_tiles_phi
577 + (iphi+_n_tiles_phi) % _n_tiles_phi;
578 }
579
580 // routines for tiled case, including some overloads of the plain
581 // BriefJet cases
582 int _tile_index(const double & eta, const double & phi) const;
583 void _tj_set_jetinfo ( TiledJet * const jet, const int _jets_index);
584 void _bj_remove_from_tiles(TiledJet * const jet);
585 void _initialise_tiles();
586 void _print_tiles(TiledJet * briefjets ) const;
587 void _add_neighbours_to_tile_union(const int tile_index,
588 std::vector<int> & tile_union, int & n_near_tiles) const;
589 void _add_untagged_neighbours_to_tile_union(const int tile_index,
590 std::vector<int> & tile_union, int & n_near_tiles);
591
592
593};
594
595
596
597//**********************************************************************
598//************** START OF INLINE MATERIAL ******************
599//**********************************************************************
600
601
602//----------------------------------------------------------------------
603// Transfer the initial jets into our internal structure
604template<class L> void ClusterSequence::_transfer_input_jets(
605 const std::vector<L> & pseudojets) {
606
607 // this will ensure that we can point to jets without difficulties
608 // arising.
609 _jets.reserve(pseudojets.size()*2);
610
611 // insert initial jets this way so that any type L that can be
612 // converted to a pseudojet will work fine (basically PseudoJet
613 // and any type that has [] subscript access to the momentum
614 // components, such as CLHEP HepLorentzVector).
615 for (unsigned int i = 0; i < pseudojets.size(); i++) {
616 _jets.push_back(pseudojets[i]);}
617
618}
619
620//----------------------------------------------------------------------
621// initialise from some generic type... Has to be made available
622// here in order for it the template aspect of it to work...
623template<class L> ClusterSequence::ClusterSequence (
624 const std::vector<L> & pseudojets,
625 const double & R,
626 const Strategy & strategy,
627 const bool & writeout_combinations) {
628
629 // transfer the initial jets (type L) into our own array
630 _transfer_input_jets(pseudojets);
631
632 // run the clustering
633 _initialise_and_run(R,strategy,writeout_combinations);
634}
635
636
637//----------------------------------------------------------------------
638/// constructor of a jet-clustering sequence from a vector of
639/// four-momenta, with the jet definition specified by jet_def
640template<class L> ClusterSequence::ClusterSequence (
641 const std::vector<L> & pseudojets,
642 const JetDefinition & jet_def,
643 const bool & writeout_combinations) {
644
645 // transfer the initial jets (type L) into our own array
646 _transfer_input_jets(pseudojets);
647
648 // run the clustering
649 _initialise_and_run(jet_def,writeout_combinations);
650}
651
652
653inline const std::vector<PseudoJet> & ClusterSequence::jets () const {
654 return _jets;
655}
656
657inline const std::vector<ClusterSequence::history_element> & ClusterSequence::history () const {
658 return _history;
659}
660
661inline unsigned int ClusterSequence::n_particles() const {return _initial_n;}
662
663
664
665//----------------------------------------------------------------------
666template <class J> inline void ClusterSequence::_bj_set_jetinfo(
667 J * const jetA, const int _jets_index) const {
668 jetA->eta = _jets[_jets_index].rap();
669 jetA->phi = _jets[_jets_index].phi_02pi();
670 jetA->kt2 = jet_scale_for_algorithm(_jets[_jets_index]);
671 jetA->_jets_index = _jets_index;
672 // initialise NN info as well
673 jetA->NN_dist = _R2;
674 jetA->NN = NULL;
675}
676
677
678
679
680//----------------------------------------------------------------------
681template <class J> inline double ClusterSequence::_bj_dist(
682 const J * const jetA, const J * const jetB) const {
683 double dphi = std::abs(jetA->phi - jetB->phi);
684 double deta = (jetA->eta - jetB->eta);
685 if (dphi > pi) {dphi = twopi - dphi;}
686 return dphi*dphi + deta*deta;
687}
688
689//----------------------------------------------------------------------
690template <class J> inline double ClusterSequence::_bj_diJ(const J * const jet) const {
691 double kt2 = jet->kt2;
692 if (jet->NN != NULL) {if (jet->NN->kt2 < kt2) {kt2 = jet->NN->kt2;}}
693 return jet->NN_dist * kt2;
694}
695
696
697//----------------------------------------------------------------------
698// set the NN for jet without checking whether in the process you might
699// have discovered a new nearest neighbour for another jet
700template <class J> inline void ClusterSequence::_bj_set_NN_nocross(
701 J * const jet, J * const head, const J * const tail) const {
702 double NN_dist = _R2;
703 J * NN = NULL;
704 if (head < jet) {
705 for (J * jetB = head; jetB != jet; jetB++) {
706 double dist = _bj_dist(jet,jetB);
707 if (dist < NN_dist) {
708 NN_dist = dist;
709 NN = jetB;
710 }
711 }
712 }
713 if (tail > jet) {
714 for (J * jetB = jet+1; jetB != tail; jetB++) {
715 double dist = _bj_dist(jet,jetB);
716 if (dist < NN_dist) {
717 NN_dist = dist;
718 NN = jetB;
719 }
720 }
721 }
722 jet->NN = NN;
723 jet->NN_dist = NN_dist;
724}
725
726
727//----------------------------------------------------------------------
728template <class J> inline void ClusterSequence::_bj_set_NN_crosscheck(J * const jet,
729 J * const head, const J * const tail) const {
730 double NN_dist = _R2;
731 J * NN = NULL;
732 for (J * jetB = head; jetB != tail; jetB++) {
733 double dist = _bj_dist(jet,jetB);
734 if (dist < NN_dist) {
735 NN_dist = dist;
736 NN = jetB;
737 }
738 if (dist < jetB->NN_dist) {
739 jetB->NN_dist = dist;
740 jetB->NN = jet;
741 }
742 }
743 jet->NN = NN;
744 jet->NN_dist = NN_dist;
745}
746
747
748
749
750FASTJET_END_NAMESPACE
751
752#endif // __FASTJET_CLUSTERSEQUENCE_HH__
Note: See TracBrowser for help on using the repository browser.