FastJet  3.1.1
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Modules Pages
ClusterSequence.cc
1 //FJSTARTHEADER
2 // $Id: ClusterSequence.cc 3685 2014-09-11 20:15:00Z salam $
3 //
4 // Copyright (c) 2005-2014, 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. They are described in the original FastJet paper,
16 // hep-ph/0512210 and in the manual, arXiv:1111.6097. If you use
17 // FastJet as part of work towards a scientific publication, please
18 // quote the version you use and include a citation to the manual and
19 // optionally also to hep-ph/0512210.
20 //
21 // FastJet is distributed in the hope that it will be useful,
22 // but WITHOUT ANY WARRANTY; without even the implied warranty of
23 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
24 // GNU General Public License for more details.
25 //
26 // You should have received a copy of the GNU General Public License
27 // along with FastJet. If not, see <http://www.gnu.org/licenses/>.
28 //----------------------------------------------------------------------
29 //FJENDHEADER
30 
31 #include "fastjet/Error.hh"
32 #include "fastjet/PseudoJet.hh"
33 #include "fastjet/ClusterSequence.hh"
34 #include "fastjet/ClusterSequenceStructure.hh"
35 #include "fastjet/version.hh" // stores the current version number
36 #include "fastjet/internal/LazyTiling9Alt.hh"
37 #include "fastjet/internal/LazyTiling9.hh"
38 #include "fastjet/internal/LazyTiling25.hh"
39 #ifndef __FJCORE__
40 #include "fastjet/internal/LazyTiling9SeparateGhosts.hh"
41 #endif // __FJCORE__
42 #include<iostream>
43 #include<sstream>
44 #include<fstream>
45 #include<cmath>
46 #include<cstdlib>
47 #include<cassert>
48 #include<string>
49 #include<set>
50 
51 FASTJET_BEGIN_NAMESPACE // defined in fastjet/internal/base.hh
52 
53 //----------------------------------------------------------------------
54 // here's where we put the main page for fastjet (as explained in the
55 // Doxygen FAQ)
56 // We put it inside the fastjet namespace to have the links without
57 // having to specify (fastjet::)
58 //......................................................................
59 /** \mainpage FastJet code documentation
60  *
61  * These pages provide automatically generated documentation for the
62  * FastJet package.
63  *
64  * \section useful_classes The most useful classes
65  *
66  * Many of the facilities of FastJet can be accessed through the three
67  * following classes:
68  *
69  * - PseudoJet: the basic class for holding the 4-momentum of a
70  * particle or a jet.
71  *
72  * - JetDefinition: the combination of a #JetAlgorithm and its
73  * associated parameters. Can also be initialised with a \ref plugins "plugin".
74  *
75  * - ClusterSequence: constructed with a vector of input (PseudoJet)
76  * particles and a JetDefinition, it computes and stores the
77  * information on how the input particles are clustered into jets.
78  *
79  * \section advanced_classes Selected more advanced classes
80  *
81  * - ClusterSequenceArea: with the help of an AreaDefinition, provides
82  * jets that also contain information about their area.
83  *
84  * \section Tools Selected additional tools
85  *
86  * - JetMedianBackgroundEstimator: with the help of a Selector, a JetDefinition and
87  * an AreaDefinition, allows one to estimate the background noise density in an event; for a simpler, quicker, effective alternative, use GridMedianBackgroundEstimator
88  *
89  * - Transformer: class from which are derived various tools for
90  * manipulating jets and accessing their substructure. Examples are
91  * Subtractor, Filter, Pruner and various taggers (e.g. JHTopTagger
92  * and MassDropTagger).
93  *
94  * \section further_info Further information
95  *
96  * - Selected classes ordered by topics can be found under the <a
97  * href="modules.html">modules</a> tab.
98  *
99  * - The complete list of classes is available under the <a
100  * href="annotated.html">classes</a> tab.
101  *
102  * - For non-class material (<a href="namespacefastjet.html#enum-members">enums</a>,
103  * <a href="namespacefastjet.html#typedef-members">typedefs</a>,
104  * <a href="namespacefastjet.html#func-members">functions</a>), see the
105  * #fastjet documentation
106  *
107  * - For further information and normal documentation, see the main <a
108  * href="http://fastjet.fr/">FastJet</a> page.
109  *
110  * \section examples Examples
111  * See our \subpage Examples page
112  */
113 
114 // define the doxygen groups
115 /// \defgroup basic_classes Fundamental FastJet classes
116 /// \defgroup area_classes Area-related classes
117 /// \defgroup sec_area_classes Secondary area-related classes
118 /// \defgroup plugins Plugins for non-native jet definitions
119 /// \defgroup selectors Selectors
120 /// \defgroup tools FastJet tools
121 /// \{ \defgroup tools_generic Generic tools
122 /// \defgroup tools_background Background subtraction
123 /// \defgroup tools_taggers Taggers
124 /// \}
125 /// \defgroup extra_info Access to extra information
126 /// \defgroup error_handling Error handling
127 /// \defgroup advanced_usage Advanced usage
128 /// \if internal_doc
129 /// \defgroup internal
130 /// \endif
131 
132 //----------------------------------------------------------------------
133 
134 
135 using namespace std;
136 
137 
138 // The following variable can be modified from within user code
139 // so as to redirect banners to an ostream other than cout.
140 //
141 // Please note that if you distribute 3rd party code
142 // that links with FastJet, that 3rd party code is NOT
143 // allowed to turn off the printing of FastJet banners
144 // by default. This requirement reflects the spirit of
145 // clause 2c of the GNU Public License (v2), under which
146 // FastJet and its plugins are distributed.
147 std::ostream * ClusterSequence::_fastjet_banner_ostr = &cout;
148 
149 
150 // destructor that guarantees proper bookkeeping for the CS Structure
151 ClusterSequence::~ClusterSequence () {
152  // set the pointer in the wrapper to this object to NULL to say that
153  // we're going out of scope
154  if (_structure_shared_ptr()){
155  ClusterSequenceStructure* csi = dynamic_cast<ClusterSequenceStructure*>(_structure_shared_ptr());
156  // normally the csi is purely internal so it really should not be
157  // NULL i.e assert should be OK
158  // (we assert rather than throw an error, since failure here is a
159  // sign of major internal problems)
160  assert(csi != NULL);
161  csi->set_associated_cs(NULL);
162 
163  // if the user had given the CS responsibility to delete itself,
164  // but then deletes the CS themselves, the following lines of
165  // code will ensure that the structure_shared_ptr will have
166  // a proper object count (so that jets associated with the CS will
167  // throw the correct error if the user tries to access their
168  // constituents).
169  if (_deletes_self_when_unused) {
170  _structure_shared_ptr.set_count(_structure_shared_ptr.use_count()
171  + _structure_use_count_after_construction);
172  }
173  }
174 }
175 
176 //-----------
177 void ClusterSequence::signal_imminent_self_deletion() const {
178  // normally if the destructor is called when
179  // _deletes_self_when_unused is true, it assumes that it's been
180  // called by the user (and it therefore resets the shared pointer
181  // count to the true count).
182  //
183  // for self deletion (called from the destructor of the CSstructure,
184  // the shared_ptr to which has just had its pointer -> 0) you do
185  // _not_ want to reset the pointer count (otherwise you will end up
186  // with a double delete on the shared pointer once you start
187  // deleting the internal structure of the CS).
188  //
189  // the following modification ensures that the count reset will not
190  // take place in the destructor
191  assert(_deletes_self_when_unused);
192  _deletes_self_when_unused = false;
193 }
194 
195 //DEP //----------------------------------------------------------------------
196 //DEP void ClusterSequence::_initialise_and_run (
197 //DEP const double R,
198 //DEP const Strategy & strategy,
199 //DEP const bool & writeout_combinations) {
200 //DEP
201 //DEP JetDefinition jet_def(_default_jet_algorithm, R, strategy);
202 //DEP _initialise_and_run(jet_def, writeout_combinations);
203 //DEP }
204 
205 
206 //----------------------------------------------------------------------
207 void ClusterSequence::_initialise_and_run (
208  const JetDefinition & jet_def_in,
209  const bool & writeout_combinations) {
210 
211  // transfer all relevant info into internal variables
212  _decant_options(jet_def_in, writeout_combinations);
213 
214  // now run
215  _initialise_and_run_no_decant();
216 }
217 
218 //----------------------------------------------------------------------
219 void ClusterSequence::_initialise_and_run_no_decant () {
220 
221  // set up the history entries for the initial particles (those
222  // currently in _jets)
223  _fill_initial_history();
224 
225  // don't run anything if the event is empty
226  if (n_particles() == 0) return;
227 
228  // ----- deal with special cases: plugins & e+e- ------
229  if (_jet_algorithm == plugin_algorithm) {
230  // allows plugin_xyz() functions to modify cluster sequence
231  _plugin_activated = true;
232  // let the plugin do its work here
233  _jet_def.plugin()->run_clustering( (*this) );
234  _plugin_activated = false;
235  _update_structure_use_count();
236  return;
237  } else if (_jet_algorithm == ee_kt_algorithm ||
238  _jet_algorithm == ee_genkt_algorithm) {
239  // ignore requested strategy
240  _strategy = N2Plain;
241  if (_jet_algorithm == ee_kt_algorithm) {
242  // make sure that R is large enough so that "beam" recomb only
243  // occurs when a single particle is left
244  // Normally, this should be automatically set to 4 from JetDefinition
245  assert(_Rparam > 2.0);
246  // this is used to renormalise the dij to get a "standard" form
247  // and our convention in e+e- will be different from that
248  // in long.inv case; NB: _invR2 name should be changed -> _renorm_dij?
249  _invR2 = 1.0;
250  } else {
251  // as of 2009-01-09, choose R to be an angular distance, in
252  // radians. Since the algorithm uses 2(1-cos(theta)) as its
253  // squared angular measure, make sure that the _R2 is defined
254  // in a similar way.
255  if (_Rparam > pi) {
256  // choose a value that ensures that back-to-back particles will
257  // always recombine
258  //_R2 = 4.0000000000001;
259  _R2 = 2 * ( 3.0 + cos(_Rparam) );
260  } else {
261  _R2 = 2 * ( 1.0 - cos(_Rparam) );
262  }
263  _invR2 = 1.0/_R2;
264  }
265  _simple_N2_cluster_EEBriefJet();
266  return;
267  } else if (_jet_algorithm == undefined_jet_algorithm) {
268  throw Error("A ClusterSequence cannot be created with an uninitialised JetDefinition");
269  }
270 
271 
272  // automatically redefine the strategy according to N if that is
273  // what the user requested -- transition points (and especially
274  // their R-dependence) are based on empirical observations for a
275  // R=0.4, 0.7 and 1.0, running on toth (3.4GHz, Pentium IV D [dual
276  // core] with 2MB of cache).
277  //-------------
278  // 2011-11-15: lowered N2Plain -> N2Tiled switchover based on some
279  // new tests on an Intel Core 2 Duo T9400 @ 2.53 GHz
280  // with 6MB cache; tests performed with lines such as
281  // ./fastjet_timing_plugins -kt -nhardest 30 -repeat 50000 -strategy -3 -R 0.5 -nev 1 < ../../data/Pythia-PtMin1000-LHC-1000ev.dat
282  if (_strategy == Best) {
283  _strategy = _best_strategy();
284 #ifdef DROP_CGAL
285  // fall back strategy for large N when CGAL is missing
286  if (_strategy == NlnN) _strategy = N2MHTLazy25;
287 #endif // DROP_CGAL
288  } else if (_strategy == BestFJ30) {
289  int N = _jets.size();
290  //if (N <= 55*max(0.5,min(1.0,_Rparam))) {// old empirical scaling with R
291  //----------------------
292  // 2011-11-15: new empirical scaling with R; NB: low-R N2Tiled
293  // could be significantly improved at low N by limiting the
294  // minimum size of tiles when R is small
295  if (min(1.0,max(0.1,_Rparam)*3.3)*N <= 30) {
296  _strategy = N2Plain;
297  } else if (N > 6200/pow(_Rparam,2.0) && _jet_def.jet_algorithm() == cambridge_algorithm) {
298  _strategy = NlnNCam;
299 #ifndef DROP_CGAL
300  } else if ((N > 16000/pow(_Rparam,1.15) && _jet_def.jet_algorithm() != antikt_algorithm)
301  || N > 35000/pow(_Rparam,1.15)) {
302  _strategy = NlnN;
303 #endif // DROP_CGAL
304  } else if (N <= 450) {
305  _strategy = N2Tiled;
306  } else {
307  _strategy = N2MinHeapTiled;
308  }
309  }
310 
311  // R >= 2pi is not supported by all clustering strategies owing to
312  // periodicity issues (a particle might cluster with itself). When
313  // R>=2pi, we therefore automatically switch to a strategy that is
314  // known to work.
315  if (_Rparam >= twopi) {
316  if ( _strategy == NlnN
317  || _strategy == NlnN3pi
318  || _strategy == NlnNCam
319  || _strategy == NlnNCam2pi2R
320  || _strategy == NlnNCam4pi) {
321 #ifdef DROP_CGAL
322  _strategy = N2MinHeapTiled;
323 #else
324  _strategy = NlnN4pi;
325 #endif
326  }
327  if (_jet_def.strategy() != Best && _strategy != _jet_def.strategy()) {
328  ostringstream oss;
329  oss << "Cluster strategy " << strategy_string(_jet_def.strategy())
330  << " automatically changed to " << strategy_string()
331  << " because the former is not supported for R = " << _Rparam
332  << " >= 2pi";
333  _changed_strategy_warning.warn(oss.str());
334  }
335  }
336 
337 
338  // run the code containing the selected strategy
339  //
340  // We order the strategies starting from the ones used by the Best
341  // strategy in the order of increasing N, then the remaining ones
342  // again in the order of increasing N.
343  if (_strategy == N2Plain) {
344  // BriefJet provides standard long.invariant kt alg.
345  this->_simple_N2_cluster_BriefJet();
346  } else if (_strategy == N2Tiled) {
347  this->_faster_tiled_N2_cluster();
348  } else if (_strategy == N2MinHeapTiled) {
349  this->_minheap_faster_tiled_N2_cluster();
350  } else if (_strategy == N2MHTLazy9Alt) {
351  // attempt to use an external tiling routine -- it manipulates
352  // the CS history via the plugin mechanism
353  _plugin_activated = true;
354  LazyTiling9Alt tiling(*this);
355  tiling.run();
356  _plugin_activated = false;
357 
358  } else if (_strategy == N2MHTLazy25) {
359  // attempt to use an external tiling routine -- it manipulates
360  // the CS history via the plugin mechanism
361  _plugin_activated = true;
362  LazyTiling25 tiling(*this);
363  tiling.run();
364  _plugin_activated = false;
365 
366  } else if (_strategy == N2MHTLazy9) {
367  // attempt to use an external tiling routine -- it manipulates
368  // the CS history via the plugin mechanism
369  _plugin_activated = true;
370  LazyTiling9 tiling(*this);
371  tiling.run();
372  _plugin_activated = false;
373 
374 #ifndef __FJCORE__
375  } else if (_strategy == N2MHTLazy9AntiKtSeparateGhosts) {
376  // attempt to use an external tiling routine -- it manipulates
377  // the CS history via the plugin mechanism
378  _plugin_activated = true;
379  LazyTiling9SeparateGhosts tiling(*this);
380  tiling.run();
381  _plugin_activated = false;
382 #else
383  throw Error("N2MHTLazy9AntiKtSeparateGhosts strategy not supported with FJCORE");
384 #endif // __FJCORE__
385 
386  } else if (_strategy == NlnN) {
387  this->_delaunay_cluster();
388  } else if (_strategy == NlnNCam) {
389  this->_CP2DChan_cluster_2piMultD();
390  } else if (_strategy == NlnN3pi || _strategy == NlnN4pi ) {
391  this->_delaunay_cluster();
392  } else if (_strategy == N3Dumb ) {
393  this->_really_dumb_cluster();
394  } else if (_strategy == N2PoorTiled) {
395  this->_tiled_N2_cluster();
396  } else if (_strategy == NlnNCam4pi) {
397  this->_CP2DChan_cluster();
398  } else if (_strategy == NlnNCam2pi2R) {
399  this->_CP2DChan_cluster_2pi2R();
400  } else {
401  ostringstream err;
402  err << "Unrecognised value for strategy: "<<_strategy;
403  throw Error(err.str());
404  }
405 
406 }
407 
408 
409 // these needs to be defined outside the class definition.
410 bool ClusterSequence::_first_time = true;
411 LimitedWarning ClusterSequence::_exclusive_warnings;
412 
413 
414 //----------------------------------------------------------------------
415 // the version string
417  return "FastJet version "+string(fastjet_version);
418 }
419 
420 
421 //----------------------------------------------------------------------
422 // prints a banner on the first call
423 void ClusterSequence::print_banner() {
424 
425  if (!_first_time) {return;}
426  _first_time = false;
427 
428  // make sure the user has not set the banner stream to NULL
429  ostream * ostr = _fastjet_banner_ostr;
430  if (!ostr) return;
431 
432  (*ostr) << "#--------------------------------------------------------------------------\n";
433  (*ostr) << "# FastJet release " << fastjet_version << endl;
434  (*ostr) << "# M. Cacciari, G.P. Salam and G. Soyez \n";
435  (*ostr) << "# A software package for jet finding and analysis at colliders \n";
436  (*ostr) << "# http://fastjet.fr \n";
437  (*ostr) << "# \n";
438  (*ostr) << "# Please cite EPJC72(2012)1896 [arXiv:1111.6097] if you use this package\n";
439  (*ostr) << "# for scientific work and optionally PLB641(2006)57 [hep-ph/0512210]. \n";
440  (*ostr) << "# \n";
441  (*ostr) << "# FastJet is provided without warranty under the terms of the GNU GPLv2.\n";
442  (*ostr) << "# It uses T. Chan's closest pair algorithm, S. Fortune's Voronoi code";
443 #ifndef DROP_CGAL
444  (*ostr) << ",\n# CGAL ";
445 #else
446  (*ostr) << "\n# ";
447 #endif // DROP_CGAL
448  (*ostr) << "and 3rd party plugin jet algorithms. See COPYING file for details.\n";
449  (*ostr) << "#--------------------------------------------------------------------------\n";
450  // make sure we really have the output done.
451  ostr->flush();
452 }
453 
454 //----------------------------------------------------------------------
455 // transfer all relevant info into internal variables
456 void ClusterSequence::_decant_options(const JetDefinition & jet_def_in,
457  const bool & writeout_combinations) {
458  // make a local copy of the jet definition (for future use)
459  _jet_def = jet_def_in;
460  _writeout_combinations = writeout_combinations;
461  // initialised the wrapper to the current CS
462  _structure_shared_ptr.reset(new ClusterSequenceStructure(this));
463 
464  _decant_options_partial();
465 }
466 
467 //----------------------------------------------------------------------
468 // transfer all relevant info into internal variables
469 void ClusterSequence::_decant_options_partial() {
470  // let the user know what's going on
471  print_banner();
472 
473  _jet_algorithm = _jet_def.jet_algorithm();
474  _Rparam = _jet_def.R(); _R2 = _Rparam*_Rparam; _invR2 = 1.0/_R2;
475  _strategy = _jet_def.strategy();
476 
477  // disallow interference from the plugin
478  _plugin_activated = false;
479 
480  // initialised the wrapper to the current CS
481  //_structure_shared_ptr.reset(new ClusterSequenceStructure(this));
482  _update_structure_use_count(); // make sure it's correct already here
483 }
484 
485 
486 //----------------------------------------------------------------------
487 // initialise the history in a standard way
488 void ClusterSequence::_fill_initial_history () {
489 
490  //if (_jets.size() == 0) {throw Error("Cannot run jet-finder on empty event");}
491 
492  // reserve sufficient space for everything
493  _jets.reserve(_jets.size()*2);
494  _history.reserve(_jets.size()*2);
495 
496  _Qtot = 0;
497 
498  for (int i = 0; i < static_cast<int>(_jets.size()) ; i++) {
499  history_element element;
500  element.parent1 = InexistentParent;
501  element.parent2 = InexistentParent;
502  element.child = Invalid;
503  element.jetp_index = i;
504  element.dij = 0.0;
505  element.max_dij_so_far = 0.0;
506 
507  _history.push_back(element);
508 
509  // do any momentum preprocessing needed by the recombination scheme
510  _jet_def.recombiner()->preprocess(_jets[i]);
511 
512  // get cross-referencing right from PseudoJets
513  _jets[i].set_cluster_hist_index(i);
514  _set_structure_shared_ptr(_jets[i]);
515 
516  // determine the total energy in the event
517  _Qtot += _jets[i].E();
518  }
519  _initial_n = _jets.size();
520  _deletes_self_when_unused = false;
521 }
522 
523 
524 //----------------------------------------------------------------------
525 string ClusterSequence::strategy_string (Strategy strategy_in) const {
526  string strategy;
527  switch(strategy_in) {
528  case NlnN:
529  strategy = "NlnN"; break;
530  case NlnN3pi:
531  strategy = "NlnN3pi"; break;
532  case NlnN4pi:
533  strategy = "NlnN4pi"; break;
534  case N2Plain:
535  strategy = "N2Plain"; break;
536  case N2Tiled:
537  strategy = "N2Tiled"; break;
538  case N2MinHeapTiled:
539  strategy = "N2MinHeapTiled"; break;
540  case N2PoorTiled:
541  strategy = "N2PoorTiled"; break;
542  case N2MHTLazy9:
543  strategy = "N2MHTLazy9"; break;
544  case N2MHTLazy9Alt:
545  strategy = "N2MHTLazy9Alt"; break;
546  case N2MHTLazy25:
547  strategy = "N2MHTLazy25"; break;
549  strategy = "N2MHTLazy9AntiKtSeparateGhosts"; break;
550  case N3Dumb:
551  strategy = "N3Dumb"; break;
552  case NlnNCam4pi:
553  strategy = "NlnNCam4pi"; break;
554  case NlnNCam2pi2R:
555  strategy = "NlnNCam2pi2R"; break;
556  case NlnNCam:
557  strategy = "NlnNCam"; break; // 2piMultD
558  case plugin_strategy:
559  strategy = "plugin strategy"; break;
560  default:
561  strategy = "Unrecognized";
562  }
563  return strategy;
564 }
565 
566 
567 double ClusterSequence::jet_scale_for_algorithm(
568  const PseudoJet & jet) const {
569  if (_jet_algorithm == kt_algorithm) {return jet.kt2();}
570  else if (_jet_algorithm == cambridge_algorithm) {return 1.0;}
571  else if (_jet_algorithm == antikt_algorithm) {
572  double kt2=jet.kt2();
573  return kt2 > 1e-300 ? 1.0/kt2 : 1e300;
574  } else if (_jet_algorithm == genkt_algorithm) {
575  double kt2 = jet.kt2();
576  double p = jet_def().extra_param();
577  if (p <= 0 && kt2 < 1e-300) kt2 = 1e-300; // dodgy safety check
578  return pow(kt2, p);
579  } else if (_jet_algorithm == cambridge_for_passive_algorithm) {
580  double kt2 = jet.kt2();
581  double lim = _jet_def.extra_param();
582  if (kt2 < lim*lim && kt2 != 0.0) {
583  return 1.0/kt2;
584  } else {return 1.0;}
585  } else {throw Error("Unrecognised jet algorithm");}
586 }
587 
588 //----------------------------------------------------------------------
589 // returns a suggestion for the best strategy to use on event
590 // multiplicity, algorithm, R, etc.
591 //
592 // Some of the work to establish the best strategy is collected in
593 // issue-tracker/2014-07-auto-strategy-selection;
594 // transition_fit_v2.fit indicates the results of the fits that we're
595 // using here. (Automatically generated by transition_fit_v2.gp).
596 //
597 // The transition to NlnN is always present, and it is the the
598 // caller's responsibility to drop back down to N2MHTLazy25 if NlnN
599 // isn't available.
600 //
601 // This routine should be called only if the jet alg is one of kt,
602 // antikt, cam or genkt.
603 Strategy ClusterSequence::_best_strategy() const {
604  int N = _jets.size();
605  // define bounded R, always above 0.1, because we don't trust any
606  // of our parametrizations below R = 0.1
607  double bounded_R = max(_Rparam, 0.1);
608 
609  // the very first test thing is a quick hard-coded test to decide
610  // if we immediately opt for N2Plain
611  if (N <= 30 || N <= 39.0/(bounded_R + 0.6)) {
612  return N2Plain;
613  }
614 
615  // Define objects that describe our various boundaries. A prefix N_
616  // indicates that boundary is for N, while L_ means it's for log(N).
617  //
618  // Hopefully having them static will ensure minimal overhead
619  // in creating them; collecting them in one place should
620  // help with updates?
621  //
622  const static _Parabola N_Tiled_to_MHT_lowR (-45.4947,54.3528,44.6283);
623  const static _Parabola L_MHT_to_MHTLazy9_lowR (0.677807,-1.05006,10.6994);
624  const static _Parabola L_MHTLazy9_to_MHTLazy25_akt_lowR(0.169967,-0.512589,12.1572);
625  const static _Parabola L_MHTLazy9_to_MHTLazy25_kt_lowR (0.16237,-0.484612,12.3373);
626  const static _Parabola L_MHTLazy9_to_MHTLazy25_cam_lowR = L_MHTLazy9_to_MHTLazy25_kt_lowR;
627  const static _Parabola L_MHTLazy25_to_NlnN_akt_lowR (0.0472051,-0.22043,15.9196);
628  const static _Parabola L_MHTLazy25_to_NlnN_kt_lowR (0.118609,-0.326811,14.8287);
629  const static _Parabola L_MHTLazy25_to_NlnN_cam_lowR (0.10119,-0.295748,14.3924);
630 
631  const static _Line L_Tiled_to_MHTLazy9_medR (-1.31304,7.29621);
632  const static _Parabola L_MHTLazy9_to_MHTLazy25_akt_medR = L_MHTLazy9_to_MHTLazy25_akt_lowR;
633  const static _Parabola L_MHTLazy9_to_MHTLazy25_kt_medR = L_MHTLazy9_to_MHTLazy25_kt_lowR;
634  const static _Parabola L_MHTLazy9_to_MHTLazy25_cam_medR = L_MHTLazy9_to_MHTLazy25_cam_lowR;
635  const static _Parabola L_MHTLazy25_to_NlnN_akt_medR = L_MHTLazy25_to_NlnN_akt_lowR;
636  const static _Parabola L_MHTLazy25_to_NlnN_kt_medR = L_MHTLazy25_to_NlnN_kt_lowR;
637  const static _Parabola L_MHTLazy25_to_NlnN_cam_medR = L_MHTLazy25_to_NlnN_cam_lowR;
638 
639  const static double N_Plain_to_MHTLazy9_largeR = 75;
640  const static double N_MHTLazy9_to_MHTLazy25_akt_largeR = 700;
641  const static double N_MHTLazy9_to_MHTLazy25_kt_largeR = 1000;
642  const static double N_MHTLazy9_to_MHTLazy25_cam_largeR = 1000;
643  const static double N_MHTLazy25_to_NlnN_akt_largeR = 100000;
644  const static double N_MHTLazy25_to_NlnN_kt_largeR = 40000;
645  const static double N_MHTLazy25_to_NlnN_cam_largeR = 15000;
646 
647  // We have timing studies only for kt, cam and antikt; for other
648  // algorithms we set the local jet_algorithm variable to the one of
649  // kt,cam,antikt that we think will be closest in behaviour to the
650  // other alg.
651  JetAlgorithm jet_algorithm;
652  if (_jet_algorithm == genkt_algorithm) {
653  // for genkt, then we set the local jet_algorithm variable (used
654  // only for strategy choice) to be either kt or antikt, depending on
655  // the p value.
656  double p = jet_def().extra_param();
657  if (p < 0.0) jet_algorithm = antikt_algorithm;
658  else jet_algorithm = kt_algorithm;
659  } else if (_jet_algorithm == cambridge_for_passive_algorithm) {
660  // we assume (but haven't tested) that using the kt-alg timing
661  // transitions should be adequate for cambridge_for_passive_algorithm
662  jet_algorithm = kt_algorithm;
663  } else {
664  jet_algorithm = _jet_algorithm;
665  }
666 
667  if (bounded_R < 0.65) {
668  // low R case
669  if (N < N_Tiled_to_MHT_lowR(bounded_R)) return N2Tiled;
670  double logN = log(double(N));
671  if (logN < L_MHT_to_MHTLazy9_lowR(bounded_R)) return N2MinHeapTiled;
672  else {
673  if (jet_algorithm == antikt_algorithm){
674  if (logN < L_MHTLazy9_to_MHTLazy25_akt_lowR(bounded_R)) return N2MHTLazy9;
675  else if (logN < L_MHTLazy25_to_NlnN_akt_lowR(bounded_R)) return N2MHTLazy25;
676  else return NlnN;
677  } else if (jet_algorithm == kt_algorithm){
678  if (logN < L_MHTLazy9_to_MHTLazy25_kt_lowR(bounded_R)) return N2MHTLazy9;
679  else if (logN < L_MHTLazy25_to_NlnN_kt_lowR(bounded_R)) return N2MHTLazy25;
680  else return NlnN;
681  } else if (jet_algorithm == cambridge_algorithm) {
682  if (logN < L_MHTLazy9_to_MHTLazy25_cam_lowR(bounded_R)) return N2MHTLazy9;
683  else if (logN < L_MHTLazy25_to_NlnN_cam_lowR(bounded_R)) return N2MHTLazy25;
684  else return NlnNCam;
685  }
686  }
687  } else if (bounded_R < 0.5*pi) {
688  // medium R case
689  double logN = log(double(N));
690  if (logN < L_Tiled_to_MHTLazy9_medR(bounded_R)) return N2Tiled;
691  else {
692  if (jet_algorithm == antikt_algorithm){
693  if (logN < L_MHTLazy9_to_MHTLazy25_akt_medR(bounded_R)) return N2MHTLazy9;
694  else if (logN < L_MHTLazy25_to_NlnN_akt_medR(bounded_R)) return N2MHTLazy25;
695  else return NlnN;
696  } else if (jet_algorithm == kt_algorithm){
697  if (logN < L_MHTLazy9_to_MHTLazy25_kt_medR(bounded_R)) return N2MHTLazy9;
698  else if (logN < L_MHTLazy25_to_NlnN_kt_medR(bounded_R)) return N2MHTLazy25;
699  else return NlnN;
700  } else if (jet_algorithm == cambridge_algorithm) {
701  if (logN < L_MHTLazy9_to_MHTLazy25_cam_medR(bounded_R)) return N2MHTLazy9;
702  else if (logN < L_MHTLazy25_to_NlnN_cam_medR(bounded_R)) return N2MHTLazy25;
703  else return NlnNCam;
704  }
705  }
706  } else {
707  // large R case (R > pi/2)
708  if (N < N_Plain_to_MHTLazy9_largeR) return N2Plain;
709  else {
710  if (jet_algorithm == antikt_algorithm){
711  if (N < N_MHTLazy9_to_MHTLazy25_akt_largeR) return N2MHTLazy9;
712  else if (N < N_MHTLazy25_to_NlnN_akt_largeR) return N2MHTLazy25;
713  else return NlnN;
714  } else if (jet_algorithm == kt_algorithm){
715  if (N < N_MHTLazy9_to_MHTLazy25_kt_largeR) return N2MHTLazy9;
716  else if (N < N_MHTLazy25_to_NlnN_kt_largeR) return N2MHTLazy25;
717  else return NlnN;
718  } else if (jet_algorithm == cambridge_algorithm) {
719  if (N < N_MHTLazy9_to_MHTLazy25_cam_largeR) return N2MHTLazy9;
720  else if (N < N_MHTLazy25_to_NlnN_cam_largeR) return N2MHTLazy25;
721  else return NlnNCam;
722  }
723  }
724  }
725 
726  bool code_should_never_reach_here = false;
727  assert(code_should_never_reach_here);
728  return N2MHTLazy9;
729 
730 }
731 
732 
733 // //----------------------------------------------------------------------
734 // /// transfer the sequence contained in other_seq into our own;
735 // /// any plugin "extras" contained in the from_seq will be lost
736 // /// from there.
737 // void ClusterSequence::transfer_from_sequence(ClusterSequence & from_seq) {
738 //
739 // if (will_delete_self_when_unused())
740 // throw(Error("cannot use CS::transfer_from_sequence after a call to delete_self_when_unused()"));
741 //
742 // // the metadata
743 // _jet_def = from_seq._jet_def ;
744 // _writeout_combinations = from_seq._writeout_combinations ;
745 // _initial_n = from_seq._initial_n ;
746 // _Rparam = from_seq._Rparam ;
747 // _R2 = from_seq._R2 ;
748 // _invR2 = from_seq._invR2 ;
749 // _strategy = from_seq._strategy ;
750 // _jet_algorithm = from_seq._jet_algorithm ;
751 // _plugin_activated = from_seq._plugin_activated ;
752 //
753 // // the data
754 // _jets = from_seq._jets;
755 // _history = from_seq._history;
756 // // the following transfers ownership of the extras from the from_seq
757 // _extras = from_seq._extras;
758 //
759 // // transfer of ownership
760 // if (_structure_shared_ptr()) {
761 // // anything that is currently associated with the cluster sequence
762 // // should be told that its cluster sequence no longer exists
763 // ClusterSequenceStructure* csi = dynamic_cast<ClusterSequenceStructure*>(_structure_shared_ptr());
764 // assert(csi != NULL);
765 // csi->set_associated_cs(NULL);
766 // }
767 // // create a new _structure_shared_ptr to reflect the fact that
768 // // this CS is essentially a new one
769 // _structure_shared_ptr.reset(new ClusterSequenceStructure(this));
770 // _update_structure_use_count();
771 //
772 // for (vector<PseudoJet>::iterator jit = _jets.begin(); jit != _jets.end(); jit++)
773 // _set_structure_shared_ptr(*jit);
774 // }
775 
776 
777 //----------------------------------------------------------------------
778 // transfer the sequence contained in other_seq into our own;
779 // any plugin "extras" contained in the from_seq will be lost
780 // from there.
781 //
782 // It also sets the ClusterSequence pointers of the PseudoJets in
783 // the history to point to this ClusterSequence
784 //
785 // The second argument is an action that will be applied on every
786 // jets in the resulting ClusterSequence
787 void ClusterSequence::transfer_from_sequence(const ClusterSequence & from_seq,
788  const FunctionOfPseudoJet<PseudoJet> * action_on_jets){
789 
790  if (will_delete_self_when_unused())
791  throw(Error("cannot use CS::transfer_from_sequence after a call to delete_self_when_unused()"));
792 
793  // the metadata
794  _jet_def = from_seq._jet_def ;
795  _writeout_combinations = from_seq._writeout_combinations ;
796  _initial_n = from_seq._initial_n ;
797  _Rparam = from_seq._Rparam ;
798  _R2 = from_seq._R2 ;
799  _invR2 = from_seq._invR2 ;
800  _strategy = from_seq._strategy ;
801  _jet_algorithm = from_seq._jet_algorithm ;
802  _plugin_activated = from_seq._plugin_activated ;
803 
804  // the data
805 
806  // apply the transformation on the jets if needed
807  if (action_on_jets)
808  _jets = (*action_on_jets)(from_seq._jets);
809  else
810  _jets = from_seq._jets;
811  _history = from_seq._history;
812  // the following shares ownership of the extras with the from_seq;
813  // no transformations will be applied to the extras
814  _extras = from_seq._extras;
815 
816  // clean up existing structure
817  if (_structure_shared_ptr()) {
818  // If there are jets associated with an old version of the CS and
819  // a new one, keeping track of when to delete the CS becomes more
820  // complex; so we don't allow this situation to occur.
821  if (_deletes_self_when_unused) throw Error("transfer_from_sequence cannot be used for a cluster sequence that deletes self when unused");
822 
823  // anything that is currently associated with the cluster sequence
824  // should be told that its cluster sequence no longer exists
825  ClusterSequenceStructure* csi = dynamic_cast<ClusterSequenceStructure*>(_structure_shared_ptr());
826  assert(csi != NULL);
827  csi->set_associated_cs(NULL);
828  }
829  // create a new _structure_shared_ptr to reflect the fact that
830  // this CS is essentially a new one
831  _structure_shared_ptr.reset(new ClusterSequenceStructure(this));
832  _update_structure_use_count();
833 
834  for (unsigned int i=0; i<_jets.size(); i++){
835  // we reset the cluster history index in case action_on_jets
836  // messed up with it
837  _jets[i].set_cluster_hist_index(from_seq._jets[i].cluster_hist_index());
838 
839  // reset the structure pointer
840  _set_structure_shared_ptr(_jets[i]);
841  }
842 }
843 
844 
845 //----------------------------------------------------------------------
846 // record an ij recombination and reset the _jets[newjet_k] momentum and
847 // user index to be those of newjet
848 void ClusterSequence::plugin_record_ij_recombination(
849  int jet_i, int jet_j, double dij,
850  const PseudoJet & newjet, int & newjet_k) {
851 
852  plugin_record_ij_recombination(jet_i, jet_j, dij, newjet_k);
853 
854  // now transfer newjet into place
855  int tmp_index = _jets[newjet_k].cluster_hist_index();
856  _jets[newjet_k] = newjet;
857  _jets[newjet_k].set_cluster_hist_index(tmp_index);
858  _set_structure_shared_ptr(_jets[newjet_k]);
859 }
860 
861 
862 //----------------------------------------------------------------------
863 // return all inclusive jets with pt > ptmin
864 vector<PseudoJet> ClusterSequence::inclusive_jets (const double ptmin) const{
865  double dcut = ptmin*ptmin;
866  int i = _history.size() - 1; // last jet
867  vector<PseudoJet> jets_local;
868  if (_jet_algorithm == kt_algorithm) {
869  while (i >= 0) {
870  // with our specific definition of dij and diB (i.e. R appears only in
871  // dij), then dij==diB is the same as the jet.perp2() and we can exploit
872  // this in selecting the jets...
873  if (_history[i].max_dij_so_far < dcut) {break;}
874  if (_history[i].parent2 == BeamJet && _history[i].dij >= dcut) {
875  // for beam jets
876  int parent1 = _history[i].parent1;
877  jets_local.push_back(_jets[_history[parent1].jetp_index]);}
878  i--;
879  }
880  } else if (_jet_algorithm == cambridge_algorithm) {
881  while (i >= 0) {
882  // inclusive jets are all at end of clustering sequence in the
883  // Cambridge algorithm -- so if we find a non-exclusive jet, then
884  // we can exit
885  if (_history[i].parent2 != BeamJet) {break;}
886  int parent1 = _history[i].parent1;
887  const PseudoJet & jet = _jets[_history[parent1].jetp_index];
888  if (jet.perp2() >= dcut) {jets_local.push_back(jet);}
889  i--;
890  }
891  } else if (_jet_algorithm == plugin_algorithm
892  || _jet_algorithm == ee_kt_algorithm
893  || _jet_algorithm == antikt_algorithm
894  || _jet_algorithm == genkt_algorithm
895  || _jet_algorithm == ee_genkt_algorithm
896  || _jet_algorithm == cambridge_for_passive_algorithm) {
897  // for inclusive jets with a plugin algorithm, we make no
898  // assumptions about anything (relation of dij to momenta,
899  // ordering of the dij, etc.)
900  while (i >= 0) {
901  if (_history[i].parent2 == BeamJet) {
902  int parent1 = _history[i].parent1;
903  const PseudoJet & jet = _jets[_history[parent1].jetp_index];
904  if (jet.perp2() >= dcut) {jets_local.push_back(jet);}
905  }
906  i--;
907  }
908  } else {throw Error("cs::inclusive_jets(...): Unrecognized jet algorithm");}
909  return jets_local;
910 }
911 
912 
913 //----------------------------------------------------------------------
914 // return the number of exclusive jets that would have been obtained
915 // running the algorithm in exclusive mode with the given dcut
916 int ClusterSequence::n_exclusive_jets (const double dcut) const {
917 
918  // first locate the point where clustering would have stopped (i.e. the
919  // first time max_dij_so_far > dcut)
920  int i = _history.size() - 1; // last jet
921  while (i >= 0) {
922  if (_history[i].max_dij_so_far <= dcut) {break;}
923  i--;
924  }
925  int stop_point = i + 1;
926  // relation between stop_point, njets assumes one extra jet disappears
927  // at each clustering.
928  int njets = 2*_initial_n - stop_point;
929  return njets;
930 }
931 
932 //----------------------------------------------------------------------
933 // return all exclusive jets that would have been obtained running
934 // the algorithm in exclusive mode with the given dcut
935 vector<PseudoJet> ClusterSequence::exclusive_jets (const double dcut) const {
936  int njets = n_exclusive_jets(dcut);
937  return exclusive_jets(njets);
938 }
939 
940 
941 //----------------------------------------------------------------------
942 // return the jets obtained by clustering the event to n jets.
943 // Throw an error if there are fewer than n particles.
944 vector<PseudoJet> ClusterSequence::exclusive_jets (const int njets) const {
945 
946  // make sure the user does not ask for more than jets than there
947  // were particles in the first place.
948  if (njets > _initial_n) {
949  ostringstream err;
950  err << "Requested " << njets << " exclusive jets, but there were only "
951  << _initial_n << " particles in the event";
952  throw Error(err.str());
953  }
954 
955  return exclusive_jets_up_to(njets);
956 }
957 
958 //----------------------------------------------------------------------
959 // return the jets obtained by clustering the event to n jets.
960 // If there are fewer than n particles, simply return all particles
961 vector<PseudoJet> ClusterSequence::exclusive_jets_up_to (const int njets) const {
962 
963  // provide a warning when extracting exclusive jets for algorithms
964  // that does not support it explicitly.
965  // Native algorithm that support it are: kt, ee_kt, Cambridge/Aachen,
966  // genkt and ee_genkt (both with p>=0)
967  // For plugins, we check Plugin::exclusive_sequence_meaningful()
968  if (( _jet_def.jet_algorithm() != kt_algorithm) &&
969  ( _jet_def.jet_algorithm() != cambridge_algorithm) &&
970  ( _jet_def.jet_algorithm() != ee_kt_algorithm) &&
971  (((_jet_def.jet_algorithm() != genkt_algorithm) &&
972  (_jet_def.jet_algorithm() != ee_genkt_algorithm)) ||
973  (_jet_def.extra_param() <0)) &&
974  ((_jet_def.jet_algorithm() != plugin_algorithm) ||
975  (!_jet_def.plugin()->exclusive_sequence_meaningful()))) {
976  _exclusive_warnings.warn("dcut and exclusive jets for jet-finders other than kt, C/A or genkt with p>=0 should be interpreted with care.");
977  }
978 
979 
980  // calculate the point where we have to stop the clustering.
981  // relation between stop_point, njets assumes one extra jet disappears
982  // at each clustering.
983  int stop_point = 2*_initial_n - njets;
984  // make sure it's safe when more jets are requested than there are particles
985  if (stop_point < _initial_n) stop_point = _initial_n;
986 
987  // some sanity checking to make sure that e+e- does not give us
988  // surprises (should we ever implement e+e-)...
989  if (2*_initial_n != static_cast<int>(_history.size())) {
990  ostringstream err;
991  err << "2*_initial_n != _history.size() -- this endangers internal assumptions!\n";
992  throw Error(err.str());
993  //assert(false);
994  }
995 
996  // now go forwards and reconstitute the jets that we have --
997  // basically for any history element, see if the parent jets to
998  // which it refers were created before the stopping point -- if they
999  // were then add them to the list, otherwise they are subsequent
1000  // recombinations of the jets that we are looking for.
1001  vector<PseudoJet> jets_local;
1002  for (unsigned int i = stop_point; i < _history.size(); i++) {
1003  int parent1 = _history[i].parent1;
1004  if (parent1 < stop_point) {
1005  jets_local.push_back(_jets[_history[parent1].jetp_index]);
1006  }
1007  int parent2 = _history[i].parent2;
1008  if (parent2 < stop_point && parent2 > 0) {
1009  jets_local.push_back(_jets[_history[parent2].jetp_index]);
1010  }
1011 
1012  }
1013 
1014  // sanity check...
1015  if (int(jets_local.size()) != min(_initial_n, njets)) {
1016  ostringstream err;
1017  err << "ClusterSequence::exclusive_jets: size of returned vector ("
1018  <<jets_local.size()<<") does not coincide with requested number of jets ("
1019  <<njets<<")";
1020  throw Error(err.str());
1021  }
1022 
1023  return jets_local;
1024 }
1025 
1026 //----------------------------------------------------------------------
1027 /// return the dmin corresponding to the recombination that went from
1028 /// n+1 to n jets
1029 double ClusterSequence::exclusive_dmerge (const int njets) const {
1030  assert(njets >= 0);
1031  if (njets >= _initial_n) {return 0.0;}
1032  return _history[2*_initial_n-njets-1].dij;
1033 }
1034 
1035 
1036 //----------------------------------------------------------------------
1037 /// return the maximum of the dmin encountered during all recombinations
1038 /// up to the one that led to an n-jet final state; identical to
1039 /// exclusive_dmerge, except in cases where the dmin do not increase
1040 /// monotonically.
1041 double ClusterSequence::exclusive_dmerge_max (const int njets) const {
1042  assert(njets >= 0);
1043  if (njets >= _initial_n) {return 0.0;}
1044  return _history[2*_initial_n-njets-1].max_dij_so_far;
1045 }
1046 
1047 
1048 //----------------------------------------------------------------------
1049 /// return a vector of all subjets of the current jet (in the sense
1050 /// of the exclusive algorithm) that would be obtained when running
1051 /// the algorithm with the given dcut.
1052 std::vector<PseudoJet> ClusterSequence::exclusive_subjets
1053  (const PseudoJet & jet, const double dcut) const {
1054 
1055  set<const history_element*> subhist;
1056 
1057  // get the set of history elements that correspond to subjets at
1058  // scale dcut
1059  get_subhist_set(subhist, jet, dcut, 0);
1060 
1061  // now transfer this into a sequence of jets
1062  vector<PseudoJet> subjets;
1063  subjets.reserve(subhist.size());
1064  for (set<const history_element*>::iterator elem = subhist.begin();
1065  elem != subhist.end(); elem++) {
1066  subjets.push_back(_jets[(*elem)->jetp_index]);
1067  }
1068  return subjets;
1069 }
1070 
1071 //----------------------------------------------------------------------
1072 /// return the size of exclusive_subjets(...); still n ln n with same
1073 /// coefficient, but marginally more efficient than manually taking
1074 /// exclusive_subjets.size()
1075 int ClusterSequence::n_exclusive_subjets(const PseudoJet & jet,
1076  const double dcut) const {
1077  set<const history_element*> subhist;
1078  // get the set of history elements that correspond to subjets at
1079  // scale dcut
1080  get_subhist_set(subhist, jet, dcut, 0);
1081  return subhist.size();
1082 }
1083 
1084 //----------------------------------------------------------------------
1085 /// return the list of subjets obtained by unclustering the supplied
1086 /// jet down to nsub subjets. Throws an error if there are fewer than
1087 /// nsub particles in the jet.
1088 std::vector<PseudoJet> ClusterSequence::exclusive_subjets
1089  (const PseudoJet & jet, int nsub) const {
1090  vector<PseudoJet> subjets = exclusive_subjets_up_to(jet, nsub);
1091  if (int(subjets.size()) < nsub) {
1092  ostringstream err;
1093  err << "Requested " << nsub << " exclusive subjets, but there were only "
1094  << subjets.size() << " particles in the jet";
1095  throw Error(err.str());
1096  }
1097  return subjets;
1098 
1099 }
1100 
1101 //----------------------------------------------------------------------
1102 /// return the list of subjets obtained by unclustering the supplied
1103 /// jet down to nsub subjets (or all constituents if there are fewer
1104 /// than nsub).
1105 std::vector<PseudoJet> ClusterSequence::exclusive_subjets_up_to
1106  (const PseudoJet & jet, int nsub) const {
1107 
1108  set<const history_element*> subhist;
1109 
1110  // prepare the vector into which we'll put the result
1111  vector<PseudoJet> subjets;
1112  if (nsub < 0) throw Error("Requested a negative number of subjets. This is nonsensical.");
1113  if (nsub == 0) return subjets;
1114 
1115  // get the set of history elements that correspond to subjets at
1116  // scale dcut
1117  get_subhist_set(subhist, jet, -1.0, nsub);
1118 
1119  // now transfer this into a sequence of jets
1120  subjets.reserve(subhist.size());
1121  for (set<const history_element*>::iterator elem = subhist.begin();
1122  elem != subhist.end(); elem++) {
1123  subjets.push_back(_jets[(*elem)->jetp_index]);
1124  }
1125  return subjets;
1126 }
1127 
1128 
1129 //----------------------------------------------------------------------
1130 /// return the dij that was present in the merging nsub+1 -> nsub
1131 /// subjets inside this jet.
1132 ///
1133 /// If the jet has nsub or fewer constituents, it will return 0.
1134 double ClusterSequence::exclusive_subdmerge(const PseudoJet & jet, int nsub) const {
1135  set<const history_element*> subhist;
1136 
1137  // get the set of history elements that correspond to subjets at
1138  // scale dcut
1139  get_subhist_set(subhist, jet, -1.0, nsub);
1140 
1141  set<const history_element*>::iterator highest = subhist.end();
1142  highest--;
1143  /// will be zero if nconst <= nsub, since highest will be an original
1144  /// particle have zero dij
1145  return (*highest)->dij;
1146 }
1147 
1148 
1149 //----------------------------------------------------------------------
1150 /// return the maximum dij that occurred in the whole event at the
1151 /// stage that the nsub+1 -> nsub merge of subjets occurred inside
1152 /// this jet.
1153 ///
1154 /// If the jet has nsub or fewer constituents, it will return 0.
1155 double ClusterSequence::exclusive_subdmerge_max(const PseudoJet & jet, int nsub) const {
1156 
1157  set<const history_element*> subhist;
1158 
1159  // get the set of history elements that correspond to subjets at
1160  // scale dcut
1161  get_subhist_set(subhist, jet, -1.0, nsub);
1162 
1163  set<const history_element*>::iterator highest = subhist.end();
1164  highest--;
1165  /// will be zero if nconst <= nsub, since highest will be an original
1166  /// particle have zero dij
1167  return (*highest)->max_dij_so_far;
1168 }
1169 
1170 
1171 
1172 //----------------------------------------------------------------------
1173 /// return a set of pointers to history entries corresponding to the
1174 /// subjets of this jet; one stops going working down through the
1175 /// subjets either when
1176 /// - there is no further to go
1177 /// - one has found maxjet entries
1178 /// - max_dij_so_far <= dcut
1179 void ClusterSequence::get_subhist_set(set<const history_element*> & subhist,
1180  const PseudoJet & jet,
1181  double dcut, int maxjet) const {
1182  assert(contains(jet));
1183 
1184  subhist.clear();
1185  subhist.insert(&(_history[jet.cluster_hist_index()]));
1186 
1187  // establish the set of jets that are relevant
1188  int njet = 1;
1189  while (true) {
1190  // first find out if we need to probe deeper into jet.
1191  // Get history element closest to end of sequence
1192  set<const history_element*>::iterator highest = subhist.end();
1193  assert (highest != subhist.begin());
1194  highest--;
1195  const history_element* elem = *highest;
1196  // make sure we haven't got too many jets
1197  if (njet == maxjet) break;
1198  // make sure it has parents
1199  if (elem->parent1 < 0) break;
1200  // make sure that we still resolve it at scale dcut
1201  if (elem->max_dij_so_far <= dcut) break;
1202 
1203  // then do so: replace "highest" with its two parents
1204  subhist.erase(highest);
1205  subhist.insert(&(_history[elem->parent1]));
1206  subhist.insert(&(_history[elem->parent2]));
1207  njet++;
1208  }
1209 }
1210 
1211 //----------------------------------------------------------------------
1212 // work through the object's history until
1213 bool ClusterSequence::object_in_jet(const PseudoJet & object,
1214  const PseudoJet & jet) const {
1215 
1216  // make sure the object conceivably belongs to this clustering
1217  // sequence
1218  assert(contains(object) && contains(jet));
1219 
1220  const PseudoJet * this_object = &object;
1221  const PseudoJet * childp;
1222  while(true) {
1223  if (this_object->cluster_hist_index() == jet.cluster_hist_index()) {
1224  return true;
1225  } else if (has_child(*this_object, childp)) {
1226  this_object = childp;
1227  } else {
1228  return false;
1229  }
1230  }
1231 }
1232 
1233 //----------------------------------------------------------------------
1234 /// if the jet has parents in the clustering, it returns true
1235 /// and sets parent1 and parent2 equal to them.
1236 ///
1237 /// if it has no parents it returns false and sets parent1 and
1238 /// parent2 to zero
1239 bool ClusterSequence::has_parents(const PseudoJet & jet, PseudoJet & parent1,
1240  PseudoJet & parent2) const {
1241 
1242  const history_element & hist = _history[jet.cluster_hist_index()];
1243 
1244  // make sure we do not run into any unexpected situations --
1245  // i.e. both parents valid, or neither
1246  assert ((hist.parent1 >= 0 && hist.parent2 >= 0) ||
1247  (hist.parent1 < 0 && hist.parent2 < 0));
1248 
1249  if (hist.parent1 < 0) {
1250  parent1 = PseudoJet(0.0,0.0,0.0,0.0);
1251  parent2 = parent1;
1252  return false;
1253  } else {
1254  parent1 = _jets[_history[hist.parent1].jetp_index];
1255  parent2 = _jets[_history[hist.parent2].jetp_index];
1256  // order the parents in decreasing pt
1257  if (parent1.perp2() < parent2.perp2()) std::swap(parent1,parent2);
1258  return true;
1259  }
1260 }
1261 
1262 //----------------------------------------------------------------------
1263 /// if the jet has a child then return true and give the child jet
1264 /// otherwise return false and set the child to zero
1265 bool ClusterSequence::has_child(const PseudoJet & jet, PseudoJet & child) const {
1266 
1267  //const history_element & hist = _history[jet.cluster_hist_index()];
1268  //
1269  //if (hist.child >= 0) {
1270  // child = _jets[_history[hist.child].jetp_index];
1271  // return true;
1272  //} else {
1273  // child = PseudoJet(0.0,0.0,0.0,0.0);
1274  // return false;
1275  //}
1276  const PseudoJet * childp;
1277  bool res = has_child(jet, childp);
1278  if (res) {
1279  child = *childp;
1280  return true;
1281  } else {
1282  child = PseudoJet(0.0,0.0,0.0,0.0);
1283  return false;
1284  }
1285 }
1286 
1287 bool ClusterSequence::has_child(const PseudoJet & jet, const PseudoJet * & childp) const {
1288 
1289  const history_element & hist = _history[jet.cluster_hist_index()];
1290 
1291  // check that this jet has a child and that the child corresponds to
1292  // a true jet [RETHINK-IF-CHANGE-NUMBERING: what is the right
1293  // behaviour if the child is the same jet but made inclusive...?]
1294  if (hist.child >= 0 && _history[hist.child].jetp_index >= 0) {
1295  childp = &(_jets[_history[hist.child].jetp_index]);
1296  return true;
1297  } else {
1298  childp = NULL;
1299  return false;
1300  }
1301 }
1302 
1303 
1304 //----------------------------------------------------------------------
1305 /// if this jet has a child (and so a partner) return true
1306 /// and give the partner, otherwise return false and set the
1307 /// partner to zero
1308 bool ClusterSequence::has_partner(const PseudoJet & jet,
1309  PseudoJet & partner) const {
1310 
1311  const history_element & hist = _history[jet.cluster_hist_index()];
1312 
1313  // make sure we have a child and that the child does not correspond
1314  // to a clustering with the beam (or some other invalid quantity)
1315  if (hist.child >= 0 && _history[hist.child].parent2 >= 0) {
1316  const history_element & child_hist = _history[hist.child];
1317  if (child_hist.parent1 == jet.cluster_hist_index()) {
1318  // partner will be child's parent2 -- for iB clustering
1319  // parent2 will not be valid
1320  partner = _jets[_history[child_hist.parent2].jetp_index];
1321  } else {
1322  // partner will be child's parent1
1323  partner = _jets[_history[child_hist.parent1].jetp_index];
1324  }
1325  return true;
1326  } else {
1327  partner = PseudoJet(0.0,0.0,0.0,0.0);
1328  return false;
1329  }
1330 }
1331 
1332 
1333 //----------------------------------------------------------------------
1334 // return a vector of the particles that make up a jet
1335 vector<PseudoJet> ClusterSequence::constituents (const PseudoJet & jet) const {
1336  vector<PseudoJet> subjets;
1337  add_constituents(jet, subjets);
1338  return subjets;
1339 }
1340 
1341 //----------------------------------------------------------------------
1342 /// output the supplied vector of jets in a format that can be read
1343 /// by an appropriate root script; the format is:
1344 /// jet-n jet-px jet-py jet-pz jet-E
1345 /// particle-n particle-rap particle-phi particle-pt
1346 /// particle-n particle-rap particle-phi particle-pt
1347 /// ...
1348 /// #END
1349 /// ... [i.e. above repeated]
1350 void ClusterSequence::print_jets_for_root(const std::vector<PseudoJet> & jets_in,
1351  ostream & ostr) const {
1352  for (unsigned i = 0; i < jets_in.size(); i++) {
1353  ostr << i << " "
1354  << jets_in[i].px() << " "
1355  << jets_in[i].py() << " "
1356  << jets_in[i].pz() << " "
1357  << jets_in[i].E() << endl;
1358  vector<PseudoJet> cst = constituents(jets_in[i]);
1359  for (unsigned j = 0; j < cst.size() ; j++) {
1360  ostr << " " << j << " "
1361  << cst[j].rap() << " "
1362  << cst[j].phi() << " "
1363  << cst[j].perp() << endl;
1364  }
1365  ostr << "#END" << endl;
1366  }
1367 }
1368 
1369 void ClusterSequence::print_jets_for_root(const std::vector<PseudoJet> & jets_in,
1370  const std::string & filename,
1371  const std::string & comment ) const {
1372  std::ofstream ostr(filename.c_str());
1373  if (comment != "") ostr << "# " << comment << endl;
1374  print_jets_for_root(jets_in, ostr);
1375 }
1376 
1377 
1378 // Not yet. Perhaps in a future release
1379 // //----------------------------------------------------------------------
1380 // // print out all inclusive jets with pt > ptmin
1381 // void ClusterSequence::print_jets (const double ptmin) const{
1382 // vector<PseudoJet> jets = sorted_by_pt(inclusive_jets(ptmin));
1383 //
1384 // for (size_t j = 0; j < jets.size(); j++) {
1385 // printf("%5u %7.3f %7.3f %9.3f\n",
1386 // j,jets[j].rap(),jets[j].phi(),jets[j].perp());
1387 // }
1388 // }
1389 
1390 //----------------------------------------------------------------------
1391 /// returns a vector of size n_particles() which indicates, for
1392 /// each of the initial particles (in the order in which they were
1393 /// supplied), which of the supplied jets it belongs to; if it does
1394 /// not belong to any of the supplied jets, the index is set to -1;
1395 vector<int> ClusterSequence::particle_jet_indices(
1396  const vector<PseudoJet> & jets_in) const {
1397 
1398  vector<int> indices(n_particles());
1399 
1400  // first label all particles as not belonging to any jets
1401  for (unsigned ipart = 0; ipart < n_particles(); ipart++)
1402  indices[ipart] = -1;
1403 
1404  // then for each of the jets relabel its consituents as belonging to
1405  // that jet
1406  for (unsigned ijet = 0; ijet < jets_in.size(); ijet++) {
1407 
1408  vector<PseudoJet> jet_constituents(constituents(jets_in[ijet]));
1409 
1410  for (unsigned ip = 0; ip < jet_constituents.size(); ip++) {
1411  // a safe (if slightly redundant) way of getting the particle
1412  // index (for initial particles it is actually safe to assume
1413  // ipart=iclust).
1414  unsigned iclust = jet_constituents[ip].cluster_hist_index();
1415  unsigned ipart = history()[iclust].jetp_index;
1416  indices[ipart] = ijet;
1417  }
1418  }
1419 
1420  return indices;
1421 }
1422 
1423 
1424 //----------------------------------------------------------------------
1425 // recursive routine that adds on constituents of jet to the subjet_vector
1426 void ClusterSequence::add_constituents (
1427  const PseudoJet & jet, vector<PseudoJet> & subjet_vector) const {
1428  // find out position in cluster history
1429  int i = jet.cluster_hist_index();
1430  int parent1 = _history[i].parent1;
1431  int parent2 = _history[i].parent2;
1432 
1433  if (parent1 == InexistentParent) {
1434  // It is an original particle (labelled by its parent having value
1435  // InexistentParent), therefore add it on to the subjet vector
1436  // Note: we add the initial particle and not simply 'jet' so that
1437  // calling add_constituents with a subtracted jet containing
1438  // only one particle will work.
1439  subjet_vector.push_back(_jets[i]);
1440  return;
1441  }
1442 
1443  // add parent 1
1444  add_constituents(_jets[_history[parent1].jetp_index], subjet_vector);
1445 
1446  // see if parent2 is a real jet; if it is then add its constituents
1447  if (parent2 != BeamJet) {
1448  add_constituents(_jets[_history[parent2].jetp_index], subjet_vector);
1449  }
1450 }
1451 
1452 
1453 
1454 //----------------------------------------------------------------------
1455 // initialise the history in a standard way
1456 void ClusterSequence::_add_step_to_history (
1457  const int step_number, const int parent1,
1458  const int parent2, const int jetp_index,
1459  const double dij) {
1460 
1461  history_element element;
1462  element.parent1 = parent1;
1463  element.parent2 = parent2;
1464  element.jetp_index = jetp_index;
1465  element.child = Invalid;
1466  element.dij = dij;
1467  element.max_dij_so_far = max(dij,_history[_history.size()-1].max_dij_so_far);
1468  _history.push_back(element);
1469 
1470  int local_step = _history.size()-1;
1471  assert(local_step == step_number);
1472 
1473  assert(parent1 >= 0);
1474  _history[parent1].child = local_step;
1475  if (parent2 >= 0) {_history[parent2].child = local_step;}
1476 
1477  // get cross-referencing right from PseudoJets
1478  if (jetp_index != Invalid) {
1479  assert(jetp_index >= 0);
1480  //cout << _jets.size() <<" "<<jetp_index<<"\n";
1481  _jets[jetp_index].set_cluster_hist_index(local_step);
1482  _set_structure_shared_ptr(_jets[jetp_index]);
1483  }
1484 
1485  if (_writeout_combinations) {
1486  cout << local_step << ": "
1487  << parent1 << " with " << parent2
1488  << "; y = "<< dij<<endl;
1489  }
1490 
1491 }
1492 
1493 
1494 
1495 
1496 //======================================================================
1497 // Return an order in which to read the history such that _history[order[i]]
1498 // will always correspond to the same set of consituent particles if
1499 // two branching histories are equivalent in terms of the particles
1500 // contained in any given pseudojet.
1501 vector<int> ClusterSequence::unique_history_order() const {
1502 
1503  // first construct an array that will tell us the lowest constituent
1504  // of a given jet -- this will always be one of the original
1505  // particles, whose order is well defined and so will help us to
1506  // follow the tree in a unique manner.
1507  valarray<int> lowest_constituent(_history.size());
1508  int hist_n = _history.size();
1509  lowest_constituent = hist_n; // give it a large number
1510  for (int i = 0; i < hist_n; i++) {
1511  // sets things up for the initial partons
1512  lowest_constituent[i] = min(lowest_constituent[i],i);
1513  // propagates them through to the children of this parton
1514  if (_history[i].child > 0) lowest_constituent[_history[i].child]
1515  = min(lowest_constituent[_history[i].child],lowest_constituent[i]);
1516  }
1517 
1518  // establish an array for what we have and have not extracted so far
1519  valarray<bool> extracted(_history.size()); extracted = false;
1520  vector<int> unique_tree;
1521  unique_tree.reserve(_history.size());
1522 
1523  // now work our way through the tree
1524  for (unsigned i = 0; i < n_particles(); i++) {
1525  if (!extracted[i]) {
1526  unique_tree.push_back(i);
1527  extracted[i] = true;
1528  _extract_tree_children(i, extracted, lowest_constituent, unique_tree);
1529  }
1530  }
1531 
1532  return unique_tree;
1533 }
1534 
1535 //======================================================================
1536 // helper for unique_history_order
1537 void ClusterSequence::_extract_tree_children(
1538  int position,
1539  valarray<bool> & extracted,
1540  const valarray<int> & lowest_constituent,
1541  vector<int> & unique_tree) const {
1542  if (!extracted[position]) {
1543  // that means we may have unidentified parents around, so go and
1544  // collect them (extracted[position]) will then be made true)
1545  _extract_tree_parents(position,extracted,lowest_constituent,unique_tree);
1546  }
1547 
1548  // now look after the children...
1549  int child = _history[position].child;
1550  if (child >= 0) _extract_tree_children(child,extracted,lowest_constituent,unique_tree);
1551 }
1552 
1553 
1554 //======================================================================
1555 // return the list of unclustered particles
1556 vector<PseudoJet> ClusterSequence::unclustered_particles() const {
1557  vector<PseudoJet> unclustered;
1558  for (unsigned i = 0; i < n_particles() ; i++) {
1559  if (_history[i].child == Invalid)
1560  unclustered.push_back(_jets[_history[i].jetp_index]);
1561  }
1562  return unclustered;
1563 }
1564 
1565 //======================================================================
1566 /// Return the list of pseudojets in the ClusterSequence that do not
1567 /// have children (and are not among the inclusive jets). They may
1568 /// result from a clustering step or may be one of the pseudojets
1569 /// returned by unclustered_particles().
1570 vector<PseudoJet> ClusterSequence::childless_pseudojets() const {
1571  vector<PseudoJet> unclustered;
1572  for (unsigned i = 0; i < _history.size() ; i++) {
1573  if ((_history[i].child == Invalid) && (_history[i].parent2 != BeamJet))
1574  unclustered.push_back(_jets[_history[i].jetp_index]);
1575  }
1576  return unclustered;
1577 }
1578 
1579 
1580 
1581 //----------------------------------------------------------------------
1582 // returns true if the cluster sequence contains this jet (i.e. jet's
1583 // structure is this cluster sequence's and the cluster history index
1584 // is in a consistent range)
1585 bool ClusterSequence::contains(const PseudoJet & jet) const {
1586  return jet.cluster_hist_index() >= 0
1587  && jet.cluster_hist_index() < int(_history.size())
1589  && jet.associated_cluster_sequence() == this;
1590 }
1591 
1592 
1593 
1594 //======================================================================
1595 // helper for unique_history_order
1596 void ClusterSequence::_extract_tree_parents(
1597  int position,
1598  valarray<bool> & extracted,
1599  const valarray<int> & lowest_constituent,
1600  vector<int> & unique_tree) const {
1601 
1602  if (!extracted[position]) {
1603  int parent1 = _history[position].parent1;
1604  int parent2 = _history[position].parent2;
1605  // where relevant order parents so that we will first treat the
1606  // one containing the smaller "lowest_constituent"
1607  if (parent1 >= 0 && parent2 >= 0) {
1608  if (lowest_constituent[parent1] > lowest_constituent[parent2])
1609  std::swap(parent1, parent2);
1610  }
1611  // then actually run through the parents to extract the constituents...
1612  if (parent1 >= 0 && !extracted[parent1])
1613  _extract_tree_parents(parent1,extracted,lowest_constituent,unique_tree);
1614  if (parent2 >= 0 && !extracted[parent2])
1615  _extract_tree_parents(parent2,extracted,lowest_constituent,unique_tree);
1616  // finally declare this position to be accounted for and push it
1617  // onto our list.
1618  unique_tree.push_back(position);
1619  extracted[position] = true;
1620  }
1621 }
1622 
1623 
1624 //======================================================================
1625 /// carries out the bookkeeping associated with the step of recombining
1626 /// jet_i and jet_j (assuming a distance dij) and returns the index
1627 /// of the recombined jet, newjet_k.
1628 void ClusterSequence::_do_ij_recombination_step(
1629  const int jet_i, const int jet_j,
1630  const double dij,
1631  int & newjet_k) {
1632 
1633  // Create the new jet by recombining the first two.
1634  //
1635  // For efficiency reasons, use a ctr that initialises only the
1636  // shared pointers, since the rest of the info will anyway be dealt
1637  // with by the recombiner.
1638  PseudoJet newjet(false);
1639  _jet_def.recombiner()->recombine(_jets[jet_i], _jets[jet_j], newjet);
1640  _jets.push_back(newjet);
1641  // original version...
1642  //_jets.push_back(_jets[jet_i] + _jets[jet_j]);
1643 
1644  // get its index
1645  newjet_k = _jets.size()-1;
1646 
1647  // get history index
1648  int newstep_k = _history.size();
1649  // and provide jet with the info
1650  _jets[newjet_k].set_cluster_hist_index(newstep_k);
1651 
1652  // finally sort out the history
1653  int hist_i = _jets[jet_i].cluster_hist_index();
1654  int hist_j = _jets[jet_j].cluster_hist_index();
1655 
1656  _add_step_to_history(newstep_k, min(hist_i, hist_j), max(hist_i,hist_j),
1657  newjet_k, dij);
1658 
1659 }
1660 
1661 
1662 //======================================================================
1663 /// carries out the bookkeeping associated with the step of recombining
1664 /// jet_i with the beam
1665 void ClusterSequence::_do_iB_recombination_step(
1666  const int jet_i, const double diB) {
1667  // get history index
1668  int newstep_k = _history.size();
1669 
1670  // recombine the jet with the beam
1671  _add_step_to_history(newstep_k,_jets[jet_i].cluster_hist_index(),BeamJet,
1672  Invalid, diB);
1673 
1674 }
1675 
1676 
1677 // make sure the static member _changed_strategy_warning is defined.
1678 LimitedWarning ClusterSequence::_changed_strategy_warning;
1679 
1680 
1681 //----------------------------------------------------------------------
1682 void ClusterSequence::_set_structure_shared_ptr(PseudoJet & j) {
1683  j.set_structure_shared_ptr(_structure_shared_ptr);
1684  // record the use count of the structure shared point to help
1685  // in case we want to ask the CS to handle its own memory
1686  _update_structure_use_count();
1687 }
1688 
1689 
1690 //----------------------------------------------------------------------
1691 void ClusterSequence::_update_structure_use_count() {
1692  // record the use count of the structure shared point to help
1693  // in case we want to ask the CS to handle its own memory
1694  _structure_use_count_after_construction = _structure_shared_ptr.use_count();
1695 }
1696 
1697 //----------------------------------------------------------------------
1698 /// by calling this routine you tell the ClusterSequence to delete
1699 /// itself when all the Pseudojets associated with it have gone out
1700 /// of scope.
1701 void ClusterSequence::delete_self_when_unused() {
1702  // the trick we use to handle this is to modify the use count;
1703  // that way the structure will be deleted when there are no external
1704  // objects left associated the CS and the structure's destructor will then
1705  // look after deleting the cluster sequence
1706 
1707  // first make sure that there is at least one other object
1708  // associated with the CS
1709  int new_count = _structure_shared_ptr.use_count() - _structure_use_count_after_construction;
1710  if (new_count <= 0) {
1711  throw Error("delete_self_when_unused may only be called if at least one object outside the CS (e.g. a jet) is already associated with the CS");
1712  }
1713 
1714  _structure_shared_ptr.set_count(new_count);
1715  _deletes_self_when_unused = true;
1716 }
1717 
1718 
1719 FASTJET_END_NAMESPACE
1720 
a version of cambridge with a special distance measure for particles whose pt is < extra_param(); thi...
like the k_t but with distance measures dij = min(kti^{2p},ktj^{2p}) Delta R_{ij}^2 / R^2 diB = 1/kti...
double kt2() const
returns the squared transverse momentum
Definition: PseudoJet.hh:145
Chan's closest pair method (in a variant with 4pi coverage), for use exclusively with the Cambridge a...
deals with clustering
best of the NlnN variants – best overall for N>10^4.
Similar to N2MHTLazy9, but uses tiles of size R/2 and a 5x5 tile grid around the particle.
legacy N ln N using 4pi coverage of cylinder
int cluster_hist_index() const
return the cluster_hist_index, intended to be used by clustering routines.
Definition: PseudoJet.hh:770
the e+e- genkt algorithm (R > 2 and p=1 gives ee_kt)
fastest from about 50..500
the longitudinally invariant kt algorithm
int jetp_index
index in _history where the current jet is recombined with another jet to form its child...
Contains any information related to the clustering that should be directly accessible to PseudoJet...
std::vector< PseudoJet > _jets
This contains the physical PseudoJets; for each PseudoJet one can find the corresponding position in ...
Chan's closest pair method (in a variant with 2pi+minimal extra variant), for use exclusively with th...
std::vector< history_element > _history
this vector will contain the branching history; for each stage, _history[i].jetp_index indicates wher...
class to provide facilities for giving warnings up to some maximum number of times and to provide glo...
any plugin algorithm supplied by the user
int parent2
index in _history where first parent of this jet was created (InexistentParent if this jet is an orig...
Like N2MHTLazy9 in a number of respects, but does not calculate ghost-ghost distances and so does not...
the plugin has been used...
virtual void set_associated_cs(const ClusterSequence *new_cs)
set the associated csw
double dij
index in the _jets vector where we will find the
void set_cluster_hist_index(const int index)
set the cluster_hist_index, intended to be used by clustering routines.
Definition: PseudoJet.hh:772
bool has_valid_cluster_sequence() const
returns true if this PseudoJet has an associated and still valid(ated) ClusterSequence.
Definition: PseudoJet.cc:447
worse even than the usual N^3 algorithms
Like to N2MHTLazy9 but uses slightly different optimizations, e.g.
base class corresponding to errors that can be thrown by FastJet
Definition: Error.hh:47
string fastjet_version_string()
return a string containing information about the release
int child
index in _history where second parent of this jet was created (InexistentParent if this jet is an ori...
automatic selection of the best (based on N), including the LazyTiled strategies that are new to FJ3...
fastest below 50
a single element in the clustering history
like the k_t but with distance measures dij = min(1/kti^2,1/ktj^2) Delta R_{ij}^2 / R^2 diB = 1/kti^2...
the e+e- kt algorithm
const ClusterSequence * associated_cluster_sequence() const
get a (const) pointer to the parent ClusterSequence (NULL if inexistent)
Definition: PseudoJet.cc:437
void set_structure_shared_ptr(const SharedPtr< PseudoJetStructureBase > &structure)
set the associated structure
Definition: PseudoJet.cc:464
the longitudinally invariant variant of the cambridge algorithm (aka Aachen algoithm).
JetAlgorithm
the various families of jet-clustering algorithm
only looks into a neighbouring tile for a particle's nearest neighbour (NN) if that particle's in-til...
Strategy
the various options for the algorithmic strategy to adopt in clustering events with kt and cambridge ...
Class to contain pseudojets, including minimal information of use to jet-clustering routines...
Definition: PseudoJet.hh:67
double max_dij_so_far
the distance corresponding to the recombination at this stage of the clustering.
legacy N ln N using 3pi coverage of cylinder.
faster that N2Tiled above about 500 particles; differs from it by retainig the di(closest j) distance...
Chan's closest pair method (in a variant with 2pi+2R coverage), for use exclusively with the Cambridg...
double perp2() const
returns the squared transverse momentum
Definition: PseudoJet.hh:141
class that is intended to hold a full definition of the jet clusterer