FastJet  3.2.1
ClusterSequence.cc
1 //FJSTARTHEADER
2 // $Id: ClusterSequence.cc 4154 2016-07-20 16:20:48Z soyez $
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.get());
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  } else if (_strategy == N2MHTLazy9AntiKtSeparateGhosts) {
375 #ifndef __FJCORE__
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 
729  assert(0 && "Code should never reach here");
730 
731  return N2MHTLazy9;
732 
733 }
734 
735 
736 // //----------------------------------------------------------------------
737 // /// transfer the sequence contained in other_seq into our own;
738 // /// any plugin "extras" contained in the from_seq will be lost
739 // /// from there.
740 // void ClusterSequence::transfer_from_sequence(ClusterSequence & from_seq) {
741 //
742 // if (will_delete_self_when_unused())
743 // throw(Error("cannot use CS::transfer_from_sequence after a call to delete_self_when_unused()"));
744 //
745 // // the metadata
746 // _jet_def = from_seq._jet_def ;
747 // _writeout_combinations = from_seq._writeout_combinations ;
748 // _initial_n = from_seq._initial_n ;
749 // _Rparam = from_seq._Rparam ;
750 // _R2 = from_seq._R2 ;
751 // _invR2 = from_seq._invR2 ;
752 // _strategy = from_seq._strategy ;
753 // _jet_algorithm = from_seq._jet_algorithm ;
754 // _plugin_activated = from_seq._plugin_activated ;
755 //
756 // // the data
757 // _jets = from_seq._jets;
758 // _history = from_seq._history;
759 // // the following transfers ownership of the extras from the from_seq
760 // _extras = from_seq._extras;
761 //
762 // // transfer of ownership
763 // if (_structure_shared_ptr()) {
764 // // anything that is currently associated with the cluster sequence
765 // // should be told that its cluster sequence no longer exists
766 // ClusterSequenceStructure* csi = dynamic_cast<ClusterSequenceStructure*>(_structure_shared_ptr());
767 // assert(csi != NULL);
768 // csi->set_associated_cs(NULL);
769 // }
770 // // create a new _structure_shared_ptr to reflect the fact that
771 // // this CS is essentially a new one
772 // _structure_shared_ptr.reset(new ClusterSequenceStructure(this));
773 // _update_structure_use_count();
774 //
775 // for (vector<PseudoJet>::iterator jit = _jets.begin(); jit != _jets.end(); jit++)
776 // _set_structure_shared_ptr(*jit);
777 // }
778 
779 
780 ClusterSequence & ClusterSequence::operator=(const ClusterSequence & cs) {
781  // self assignment is trivial
782  if (&cs != this) {
783  _deletes_self_when_unused = false;
784  transfer_from_sequence(cs);
785  }
786  return *this;
787 }
788 
789 //----------------------------------------------------------------------
790 // transfer the sequence contained in other_seq into our own;
791 // any plugin "extras" contained in the from_seq will be lost
792 // from there.
793 //
794 // It also sets the ClusterSequence pointers of the PseudoJets in
795 // the history to point to this ClusterSequence
796 //
797 // The second argument is an action that will be applied on every
798 // jets in the resulting ClusterSequence
799 void ClusterSequence::transfer_from_sequence(const ClusterSequence & from_seq,
800  const FunctionOfPseudoJet<PseudoJet> * action_on_jets){
801 
802  if (will_delete_self_when_unused())
803  throw(Error("cannot use CS::transfer_from_sequence after a call to delete_self_when_unused()"));
804 
805  // the metadata
806  _jet_def = from_seq._jet_def ;
807  _writeout_combinations = from_seq._writeout_combinations ;
808  _initial_n = from_seq._initial_n ;
809  _Rparam = from_seq._Rparam ;
810  _R2 = from_seq._R2 ;
811  _invR2 = from_seq._invR2 ;
812  _strategy = from_seq._strategy ;
813  _jet_algorithm = from_seq._jet_algorithm ;
814  _plugin_activated = from_seq._plugin_activated ;
815 
816  // the data
817 
818  // apply the transformation on the jets if needed
819  if (action_on_jets)
820  _jets = (*action_on_jets)(from_seq._jets);
821  else
822  _jets = from_seq._jets;
823  _history = from_seq._history;
824  // the following shares ownership of the extras with the from_seq;
825  // no transformations will be applied to the extras
826  _extras = from_seq._extras;
827 
828  // clean up existing structure
829  if (_structure_shared_ptr) {
830  // If there are jets associated with an old version of the CS and
831  // a new one, keeping track of when to delete the CS becomes more
832  // complex; so we don't allow this situation to occur.
833  if (_deletes_self_when_unused) throw Error("transfer_from_sequence cannot be used for a cluster sequence that deletes self when unused");
834 
835  // anything that is currently associated with the cluster sequence
836  // should be told that its cluster sequence no longer exists
837  ClusterSequenceStructure* csi = dynamic_cast<ClusterSequenceStructure*>(_structure_shared_ptr.get());
838  assert(csi != NULL);
839  csi->set_associated_cs(NULL);
840  }
841  // create a new _structure_shared_ptr to reflect the fact that
842  // this CS is essentially a new one
843  _structure_shared_ptr.reset(new ClusterSequenceStructure(this));
844  _update_structure_use_count();
845 
846  for (unsigned int i=0; i<_jets.size(); i++){
847  // we reset the cluster history index in case action_on_jets
848  // messed up with it
849  _jets[i].set_cluster_hist_index(from_seq._jets[i].cluster_hist_index());
850 
851  // reset the structure pointer
852  _set_structure_shared_ptr(_jets[i]);
853  }
854 }
855 
856 
857 //----------------------------------------------------------------------
858 // record an ij recombination and reset the _jets[newjet_k] momentum and
859 // user index to be those of newjet
860 void ClusterSequence::plugin_record_ij_recombination(
861  int jet_i, int jet_j, double dij,
862  const PseudoJet & newjet, int & newjet_k) {
863 
864  plugin_record_ij_recombination(jet_i, jet_j, dij, newjet_k);
865 
866  // now transfer newjet into place
867  int tmp_index = _jets[newjet_k].cluster_hist_index();
868  _jets[newjet_k] = newjet;
869  _jets[newjet_k].set_cluster_hist_index(tmp_index);
870  _set_structure_shared_ptr(_jets[newjet_k]);
871 }
872 
873 
874 //----------------------------------------------------------------------
875 // return all inclusive jets with pt > ptmin
876 vector<PseudoJet> ClusterSequence::inclusive_jets (const double ptmin) const{
877  double dcut = ptmin*ptmin;
878  int i = _history.size() - 1; // last jet
879  vector<PseudoJet> jets_local;
880  if (_jet_algorithm == kt_algorithm) {
881  while (i >= 0) {
882  // with our specific definition of dij and diB (i.e. R appears only in
883  // dij), then dij==diB is the same as the jet.perp2() and we can exploit
884  // this in selecting the jets...
885  if (_history[i].max_dij_so_far < dcut) {break;}
886  if (_history[i].parent2 == BeamJet && _history[i].dij >= dcut) {
887  // for beam jets
888  int parent1 = _history[i].parent1;
889  jets_local.push_back(_jets[_history[parent1].jetp_index]);}
890  i--;
891  }
892  } else if (_jet_algorithm == cambridge_algorithm) {
893  while (i >= 0) {
894  // inclusive jets are all at end of clustering sequence in the
895  // Cambridge algorithm -- so if we find a non-exclusive jet, then
896  // we can exit
897  if (_history[i].parent2 != BeamJet) {break;}
898  int parent1 = _history[i].parent1;
899  const PseudoJet & jet = _jets[_history[parent1].jetp_index];
900  if (jet.perp2() >= dcut) {jets_local.push_back(jet);}
901  i--;
902  }
903  } else if (_jet_algorithm == plugin_algorithm
904  || _jet_algorithm == ee_kt_algorithm
905  || _jet_algorithm == antikt_algorithm
906  || _jet_algorithm == genkt_algorithm
907  || _jet_algorithm == ee_genkt_algorithm
908  || _jet_algorithm == cambridge_for_passive_algorithm) {
909  // for inclusive jets with a plugin algorithm, we make no
910  // assumptions about anything (relation of dij to momenta,
911  // ordering of the dij, etc.)
912  while (i >= 0) {
913  if (_history[i].parent2 == BeamJet) {
914  int parent1 = _history[i].parent1;
915  const PseudoJet & jet = _jets[_history[parent1].jetp_index];
916  if (jet.perp2() >= dcut) {jets_local.push_back(jet);}
917  }
918  i--;
919  }
920  } else {throw Error("cs::inclusive_jets(...): Unrecognized jet algorithm");}
921  return jets_local;
922 }
923 
924 
925 //----------------------------------------------------------------------
926 // return the number of exclusive jets that would have been obtained
927 // running the algorithm in exclusive mode with the given dcut
928 int ClusterSequence::n_exclusive_jets (const double dcut) const {
929 
930  // first locate the point where clustering would have stopped (i.e. the
931  // first time max_dij_so_far > dcut)
932  int i = _history.size() - 1; // last jet
933  while (i >= 0) {
934  if (_history[i].max_dij_so_far <= dcut) {break;}
935  i--;
936  }
937  int stop_point = i + 1;
938  // relation between stop_point, njets assumes one extra jet disappears
939  // at each clustering.
940  int njets = 2*_initial_n - stop_point;
941  return njets;
942 }
943 
944 //----------------------------------------------------------------------
945 // return all exclusive jets that would have been obtained running
946 // the algorithm in exclusive mode with the given dcut
947 vector<PseudoJet> ClusterSequence::exclusive_jets (const double dcut) const {
948  int njets = n_exclusive_jets(dcut);
949  return exclusive_jets(njets);
950 }
951 
952 
953 //----------------------------------------------------------------------
954 // return the jets obtained by clustering the event to n jets.
955 // Throw an error if there are fewer than n particles.
956 vector<PseudoJet> ClusterSequence::exclusive_jets (const int njets) const {
957 
958  // make sure the user does not ask for more than jets than there
959  // were particles in the first place.
960  if (njets > _initial_n) {
961  ostringstream err;
962  err << "Requested " << njets << " exclusive jets, but there were only "
963  << _initial_n << " particles in the event";
964  throw Error(err.str());
965  }
966 
967  return exclusive_jets_up_to(njets);
968 }
969 
970 //----------------------------------------------------------------------
971 // return the jets obtained by clustering the event to n jets.
972 // If there are fewer than n particles, simply return all particles
973 vector<PseudoJet> ClusterSequence::exclusive_jets_up_to (const int njets) const {
974 
975  // provide a warning when extracting exclusive jets for algorithms
976  // that does not support it explicitly.
977  // Native algorithm that support it are: kt, ee_kt, Cambridge/Aachen,
978  // genkt and ee_genkt (both with p>=0)
979  // For plugins, we check Plugin::exclusive_sequence_meaningful()
980  if (( _jet_def.jet_algorithm() != kt_algorithm) &&
981  ( _jet_def.jet_algorithm() != cambridge_algorithm) &&
982  ( _jet_def.jet_algorithm() != ee_kt_algorithm) &&
983  (((_jet_def.jet_algorithm() != genkt_algorithm) &&
984  (_jet_def.jet_algorithm() != ee_genkt_algorithm)) ||
985  (_jet_def.extra_param() <0)) &&
986  ((_jet_def.jet_algorithm() != plugin_algorithm) ||
987  (!_jet_def.plugin()->exclusive_sequence_meaningful()))) {
988  _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.");
989  }
990 
991 
992  // calculate the point where we have to stop the clustering.
993  // relation between stop_point, njets assumes one extra jet disappears
994  // at each clustering.
995  int stop_point = 2*_initial_n - njets;
996  // make sure it's safe when more jets are requested than there are particles
997  if (stop_point < _initial_n) stop_point = _initial_n;
998 
999  // some sanity checking to make sure that e+e- does not give us
1000  // surprises (should we ever implement e+e-)...
1001  if (2*_initial_n != static_cast<int>(_history.size())) {
1002  ostringstream err;
1003  err << "2*_initial_n != _history.size() -- this endangers internal assumptions!\n";
1004  throw Error(err.str());
1005  //assert(false);
1006  }
1007 
1008  // now go forwards and reconstitute the jets that we have --
1009  // basically for any history element, see if the parent jets to
1010  // which it refers were created before the stopping point -- if they
1011  // were then add them to the list, otherwise they are subsequent
1012  // recombinations of the jets that we are looking for.
1013  vector<PseudoJet> jets_local;
1014  for (unsigned int i = stop_point; i < _history.size(); i++) {
1015  int parent1 = _history[i].parent1;
1016  if (parent1 < stop_point) {
1017  jets_local.push_back(_jets[_history[parent1].jetp_index]);
1018  }
1019  int parent2 = _history[i].parent2;
1020  if (parent2 < stop_point && parent2 > 0) {
1021  jets_local.push_back(_jets[_history[parent2].jetp_index]);
1022  }
1023 
1024  }
1025 
1026  // sanity check...
1027  if (int(jets_local.size()) != min(_initial_n, njets)) {
1028  ostringstream err;
1029  err << "ClusterSequence::exclusive_jets: size of returned vector ("
1030  <<jets_local.size()<<") does not coincide with requested number of jets ("
1031  <<njets<<")";
1032  throw Error(err.str());
1033  }
1034 
1035  return jets_local;
1036 }
1037 
1038 //----------------------------------------------------------------------
1039 /// return the dmin corresponding to the recombination that went from
1040 /// n+1 to n jets
1041 double ClusterSequence::exclusive_dmerge (const int njets) const {
1042  assert(njets >= 0);
1043  if (njets >= _initial_n) {return 0.0;}
1044  return _history[2*_initial_n-njets-1].dij;
1045 }
1046 
1047 
1048 //----------------------------------------------------------------------
1049 /// return the maximum of the dmin encountered during all recombinations
1050 /// up to the one that led to an n-jet final state; identical to
1051 /// exclusive_dmerge, except in cases where the dmin do not increase
1052 /// monotonically.
1053 double ClusterSequence::exclusive_dmerge_max (const int njets) const {
1054  assert(njets >= 0);
1055  if (njets >= _initial_n) {return 0.0;}
1056  return _history[2*_initial_n-njets-1].max_dij_so_far;
1057 }
1058 
1059 
1060 //----------------------------------------------------------------------
1061 /// return a vector of all subjets of the current jet (in the sense
1062 /// of the exclusive algorithm) that would be obtained when running
1063 /// the algorithm with the given dcut.
1064 std::vector<PseudoJet> ClusterSequence::exclusive_subjets
1065  (const PseudoJet & jet, const double dcut) const {
1066 
1067  set<const history_element*> subhist;
1068 
1069  // get the set of history elements that correspond to subjets at
1070  // scale dcut
1071  get_subhist_set(subhist, jet, dcut, 0);
1072 
1073  // now transfer this into a sequence of jets
1074  vector<PseudoJet> subjets;
1075  subjets.reserve(subhist.size());
1076  for (set<const history_element*>::iterator elem = subhist.begin();
1077  elem != subhist.end(); elem++) {
1078  subjets.push_back(_jets[(*elem)->jetp_index]);
1079  }
1080  return subjets;
1081 }
1082 
1083 //----------------------------------------------------------------------
1084 /// return the size of exclusive_subjets(...); still n ln n with same
1085 /// coefficient, but marginally more efficient than manually taking
1086 /// exclusive_subjets.size()
1087 int ClusterSequence::n_exclusive_subjets(const PseudoJet & jet,
1088  const double dcut) const {
1089  set<const history_element*> subhist;
1090  // get the set of history elements that correspond to subjets at
1091  // scale dcut
1092  get_subhist_set(subhist, jet, dcut, 0);
1093  return subhist.size();
1094 }
1095 
1096 //----------------------------------------------------------------------
1097 /// return the list of subjets obtained by unclustering the supplied
1098 /// jet down to nsub subjets. Throws an error if there are fewer than
1099 /// nsub particles in the jet.
1100 std::vector<PseudoJet> ClusterSequence::exclusive_subjets
1101  (const PseudoJet & jet, int nsub) const {
1102  vector<PseudoJet> subjets = exclusive_subjets_up_to(jet, nsub);
1103  if (int(subjets.size()) < nsub) {
1104  ostringstream err;
1105  err << "Requested " << nsub << " exclusive subjets, but there were only "
1106  << subjets.size() << " particles in the jet";
1107  throw Error(err.str());
1108  }
1109  return subjets;
1110 
1111 }
1112 
1113 //----------------------------------------------------------------------
1114 /// return the list of subjets obtained by unclustering the supplied
1115 /// jet down to nsub subjets (or all constituents if there are fewer
1116 /// than nsub).
1117 std::vector<PseudoJet> ClusterSequence::exclusive_subjets_up_to
1118  (const PseudoJet & jet, int nsub) const {
1119 
1120  set<const history_element*> subhist;
1121 
1122  // prepare the vector into which we'll put the result
1123  vector<PseudoJet> subjets;
1124  if (nsub < 0) throw Error("Requested a negative number of subjets. This is nonsensical.");
1125  if (nsub == 0) return subjets;
1126 
1127  // get the set of history elements that correspond to subjets at
1128  // scale dcut
1129  get_subhist_set(subhist, jet, -1.0, nsub);
1130 
1131  // now transfer this into a sequence of jets
1132  subjets.reserve(subhist.size());
1133  for (set<const history_element*>::iterator elem = subhist.begin();
1134  elem != subhist.end(); elem++) {
1135  subjets.push_back(_jets[(*elem)->jetp_index]);
1136  }
1137  return subjets;
1138 }
1139 
1140 
1141 //----------------------------------------------------------------------
1142 /// return the dij that was present in the merging nsub+1 -> nsub
1143 /// subjets inside this jet.
1144 ///
1145 /// If the jet has nsub or fewer constituents, it will return 0.
1146 double ClusterSequence::exclusive_subdmerge(const PseudoJet & jet, int nsub) const {
1147  set<const history_element*> subhist;
1148 
1149  // get the set of history elements that correspond to subjets at
1150  // scale dcut
1151  get_subhist_set(subhist, jet, -1.0, nsub);
1152 
1153  set<const history_element*>::iterator highest = subhist.end();
1154  highest--;
1155  /// will be zero if nconst <= nsub, since highest will be an original
1156  /// particle have zero dij
1157  return (*highest)->dij;
1158 }
1159 
1160 
1161 //----------------------------------------------------------------------
1162 /// return the maximum dij that occurred in the whole event at the
1163 /// stage that the nsub+1 -> nsub merge of subjets occurred inside
1164 /// this jet.
1165 ///
1166 /// If the jet has nsub or fewer constituents, it will return 0.
1167 double ClusterSequence::exclusive_subdmerge_max(const PseudoJet & jet, int nsub) const {
1168 
1169  set<const history_element*> subhist;
1170 
1171  // get the set of history elements that correspond to subjets at
1172  // scale dcut
1173  get_subhist_set(subhist, jet, -1.0, nsub);
1174 
1175  set<const history_element*>::iterator highest = subhist.end();
1176  highest--;
1177  /// will be zero if nconst <= nsub, since highest will be an original
1178  /// particle have zero dij
1179  return (*highest)->max_dij_so_far;
1180 }
1181 
1182 
1183 
1184 //----------------------------------------------------------------------
1185 /// return a set of pointers to history entries corresponding to the
1186 /// subjets of this jet; one stops going working down through the
1187 /// subjets either when
1188 /// - there is no further to go
1189 /// - one has found maxjet entries
1190 /// - max_dij_so_far <= dcut
1191 void ClusterSequence::get_subhist_set(set<const history_element*> & subhist,
1192  const PseudoJet & jet,
1193  double dcut, int maxjet) const {
1194  assert(contains(jet));
1195 
1196  subhist.clear();
1197  subhist.insert(&(_history[jet.cluster_hist_index()]));
1198 
1199  // establish the set of jets that are relevant
1200  int njet = 1;
1201  while (true) {
1202  // first find out if we need to probe deeper into jet.
1203  // Get history element closest to end of sequence
1204  set<const history_element*>::iterator highest = subhist.end();
1205  assert (highest != subhist.begin());
1206  highest--;
1207  const history_element* elem = *highest;
1208  // make sure we haven't got too many jets
1209  if (njet == maxjet) break;
1210  // make sure it has parents
1211  if (elem->parent1 < 0) break;
1212  // make sure that we still resolve it at scale dcut
1213  if (elem->max_dij_so_far <= dcut) break;
1214 
1215  // then do so: replace "highest" with its two parents
1216  subhist.erase(highest);
1217  subhist.insert(&(_history[elem->parent1]));
1218  subhist.insert(&(_history[elem->parent2]));
1219  njet++;
1220  }
1221 }
1222 
1223 //----------------------------------------------------------------------
1224 // work through the object's history until
1225 bool ClusterSequence::object_in_jet(const PseudoJet & object,
1226  const PseudoJet & jet) const {
1227 
1228  // make sure the object conceivably belongs to this clustering
1229  // sequence
1230  assert(contains(object) && contains(jet));
1231 
1232  const PseudoJet * this_object = &object;
1233  const PseudoJet * childp;
1234  while(true) {
1235  if (this_object->cluster_hist_index() == jet.cluster_hist_index()) {
1236  return true;
1237  } else if (has_child(*this_object, childp)) {
1238  this_object = childp;
1239  } else {
1240  return false;
1241  }
1242  }
1243 }
1244 
1245 //----------------------------------------------------------------------
1246 /// if the jet has parents in the clustering, it returns true
1247 /// and sets parent1 and parent2 equal to them.
1248 ///
1249 /// if it has no parents it returns false and sets parent1 and
1250 /// parent2 to zero
1251 bool ClusterSequence::has_parents(const PseudoJet & jet, PseudoJet & parent1,
1252  PseudoJet & parent2) const {
1253 
1254  const history_element & hist = _history[jet.cluster_hist_index()];
1255 
1256  // make sure we do not run into any unexpected situations --
1257  // i.e. both parents valid, or neither
1258  assert ((hist.parent1 >= 0 && hist.parent2 >= 0) ||
1259  (hist.parent1 < 0 && hist.parent2 < 0));
1260 
1261  if (hist.parent1 < 0) {
1262  parent1 = PseudoJet(0.0,0.0,0.0,0.0);
1263  parent2 = parent1;
1264  return false;
1265  } else {
1266  parent1 = _jets[_history[hist.parent1].jetp_index];
1267  parent2 = _jets[_history[hist.parent2].jetp_index];
1268  // order the parents in decreasing pt
1269  if (parent1.perp2() < parent2.perp2()) std::swap(parent1,parent2);
1270  return true;
1271  }
1272 }
1273 
1274 //----------------------------------------------------------------------
1275 /// if the jet has a child then return true and give the child jet
1276 /// otherwise return false and set the child to zero
1277 bool ClusterSequence::has_child(const PseudoJet & jet, PseudoJet & child) const {
1278 
1279  //const history_element & hist = _history[jet.cluster_hist_index()];
1280  //
1281  //if (hist.child >= 0) {
1282  // child = _jets[_history[hist.child].jetp_index];
1283  // return true;
1284  //} else {
1285  // child = PseudoJet(0.0,0.0,0.0,0.0);
1286  // return false;
1287  //}
1288  const PseudoJet * childp;
1289  bool res = has_child(jet, childp);
1290  if (res) {
1291  child = *childp;
1292  return true;
1293  } else {
1294  child = PseudoJet(0.0,0.0,0.0,0.0);
1295  return false;
1296  }
1297 }
1298 
1299 bool ClusterSequence::has_child(const PseudoJet & jet, const PseudoJet * & childp) const {
1300 
1301  const history_element & hist = _history[jet.cluster_hist_index()];
1302 
1303  // check that this jet has a child and that the child corresponds to
1304  // a true jet [RETHINK-IF-CHANGE-NUMBERING: what is the right
1305  // behaviour if the child is the same jet but made inclusive...?]
1306  if (hist.child >= 0 && _history[hist.child].jetp_index >= 0) {
1307  childp = &(_jets[_history[hist.child].jetp_index]);
1308  return true;
1309  } else {
1310  childp = NULL;
1311  return false;
1312  }
1313 }
1314 
1315 
1316 //----------------------------------------------------------------------
1317 /// if this jet has a child (and so a partner) return true
1318 /// and give the partner, otherwise return false and set the
1319 /// partner to zero
1320 bool ClusterSequence::has_partner(const PseudoJet & jet,
1321  PseudoJet & partner) const {
1322 
1323  const history_element & hist = _history[jet.cluster_hist_index()];
1324 
1325  // make sure we have a child and that the child does not correspond
1326  // to a clustering with the beam (or some other invalid quantity)
1327  if (hist.child >= 0 && _history[hist.child].parent2 >= 0) {
1328  const history_element & child_hist = _history[hist.child];
1329  if (child_hist.parent1 == jet.cluster_hist_index()) {
1330  // partner will be child's parent2 -- for iB clustering
1331  // parent2 will not be valid
1332  partner = _jets[_history[child_hist.parent2].jetp_index];
1333  } else {
1334  // partner will be child's parent1
1335  partner = _jets[_history[child_hist.parent1].jetp_index];
1336  }
1337  return true;
1338  } else {
1339  partner = PseudoJet(0.0,0.0,0.0,0.0);
1340  return false;
1341  }
1342 }
1343 
1344 
1345 //----------------------------------------------------------------------
1346 // return a vector of the particles that make up a jet
1347 vector<PseudoJet> ClusterSequence::constituents (const PseudoJet & jet) const {
1348  vector<PseudoJet> subjets;
1349  add_constituents(jet, subjets);
1350  return subjets;
1351 }
1352 
1353 //----------------------------------------------------------------------
1354 /// output the supplied vector of jets in a format that can be read
1355 /// by an appropriate root script; the format is:
1356 /// jet-n jet-px jet-py jet-pz jet-E
1357 /// particle-n particle-rap particle-phi particle-pt
1358 /// particle-n particle-rap particle-phi particle-pt
1359 /// ...
1360 /// #END
1361 /// ... [i.e. above repeated]
1362 void ClusterSequence::print_jets_for_root(const std::vector<PseudoJet> & jets_in,
1363  ostream & ostr) const {
1364  for (unsigned i = 0; i < jets_in.size(); i++) {
1365  ostr << i << " "
1366  << jets_in[i].px() << " "
1367  << jets_in[i].py() << " "
1368  << jets_in[i].pz() << " "
1369  << jets_in[i].E() << endl;
1370  vector<PseudoJet> cst = constituents(jets_in[i]);
1371  for (unsigned j = 0; j < cst.size() ; j++) {
1372  ostr << " " << j << " "
1373  << cst[j].rap() << " "
1374  << cst[j].phi() << " "
1375  << cst[j].perp() << endl;
1376  }
1377  ostr << "#END" << endl;
1378  }
1379 }
1380 
1381 void ClusterSequence::print_jets_for_root(const std::vector<PseudoJet> & jets_in,
1382  const std::string & filename,
1383  const std::string & comment ) const {
1384  std::ofstream ostr(filename.c_str());
1385  if (comment != "") ostr << "# " << comment << endl;
1386  print_jets_for_root(jets_in, ostr);
1387 }
1388 
1389 
1390 // Not yet. Perhaps in a future release
1391 // //----------------------------------------------------------------------
1392 // // print out all inclusive jets with pt > ptmin
1393 // void ClusterSequence::print_jets (const double ptmin) const{
1394 // vector<PseudoJet> jets = sorted_by_pt(inclusive_jets(ptmin));
1395 //
1396 // for (size_t j = 0; j < jets.size(); j++) {
1397 // printf("%5u %7.3f %7.3f %9.3f\n",
1398 // j,jets[j].rap(),jets[j].phi(),jets[j].perp());
1399 // }
1400 // }
1401 
1402 //----------------------------------------------------------------------
1403 /// returns a vector of size n_particles() which indicates, for
1404 /// each of the initial particles (in the order in which they were
1405 /// supplied), which of the supplied jets it belongs to; if it does
1406 /// not belong to any of the supplied jets, the index is set to -1;
1407 vector<int> ClusterSequence::particle_jet_indices(
1408  const vector<PseudoJet> & jets_in) const {
1409 
1410  vector<int> indices(n_particles());
1411 
1412  // first label all particles as not belonging to any jets
1413  for (unsigned ipart = 0; ipart < n_particles(); ipart++)
1414  indices[ipart] = -1;
1415 
1416  // then for each of the jets relabel its consituents as belonging to
1417  // that jet
1418  for (unsigned ijet = 0; ijet < jets_in.size(); ijet++) {
1419 
1420  vector<PseudoJet> jet_constituents(constituents(jets_in[ijet]));
1421 
1422  for (unsigned ip = 0; ip < jet_constituents.size(); ip++) {
1423  // a safe (if slightly redundant) way of getting the particle
1424  // index (for initial particles it is actually safe to assume
1425  // ipart=iclust).
1426  unsigned iclust = jet_constituents[ip].cluster_hist_index();
1427  unsigned ipart = history()[iclust].jetp_index;
1428  indices[ipart] = ijet;
1429  }
1430  }
1431 
1432  return indices;
1433 }
1434 
1435 
1436 //----------------------------------------------------------------------
1437 // recursive routine that adds on constituents of jet to the subjet_vector
1438 void ClusterSequence::add_constituents (
1439  const PseudoJet & jet, vector<PseudoJet> & subjet_vector) const {
1440  // find out position in cluster history
1441  int i = jet.cluster_hist_index();
1442  int parent1 = _history[i].parent1;
1443  int parent2 = _history[i].parent2;
1444 
1445  if (parent1 == InexistentParent) {
1446  // It is an original particle (labelled by its parent having value
1447  // InexistentParent), therefore add it on to the subjet vector
1448  // Note: we add the initial particle and not simply 'jet' so that
1449  // calling add_constituents with a subtracted jet containing
1450  // only one particle will work.
1451  subjet_vector.push_back(_jets[i]);
1452  return;
1453  }
1454 
1455  // add parent 1
1456  add_constituents(_jets[_history[parent1].jetp_index], subjet_vector);
1457 
1458  // see if parent2 is a real jet; if it is then add its constituents
1459  if (parent2 != BeamJet) {
1460  add_constituents(_jets[_history[parent2].jetp_index], subjet_vector);
1461  }
1462 }
1463 
1464 
1465 
1466 //----------------------------------------------------------------------
1467 // initialise the history in a standard way
1468 void ClusterSequence::_add_step_to_history (
1469  //NO_LONGER_USED: const int step_number,
1470  const int parent1,
1471  const int parent2, const int jetp_index,
1472  const double dij) {
1473 
1474  history_element element;
1475  element.parent1 = parent1;
1476  element.parent2 = parent2;
1477  element.jetp_index = jetp_index;
1478  element.child = Invalid;
1479  element.dij = dij;
1480  element.max_dij_so_far = max(dij,_history[_history.size()-1].max_dij_so_far);
1481  _history.push_back(element);
1482 
1483  int local_step = _history.size()-1;
1484  //#ifndef __NO_ASSERTS__
1485  //assert(local_step == step_number);
1486  //#endif
1487 
1488  // sanity check: make sure the particles have not already been recombined
1489  //
1490  // Note that good practice would make this an assert (since this is
1491  // a serious internal issue). However, we decided to throw an
1492  // InternalError so that the end user can decide to catch it and
1493  // retry the clustering with a different strategy.
1494 
1495  assert(parent1 >= 0);
1496  if (_history[parent1].child != Invalid){
1497  throw InternalError("trying to recomine an object that has previsously been recombined");
1498  }
1499  _history[parent1].child = local_step;
1500  if (parent2 >= 0) {
1501  if (_history[parent2].child != Invalid){
1502  throw InternalError("trying to recomine an object that has previsously been recombined");
1503  }
1504  _history[parent2].child = local_step;
1505  }
1506 
1507  // get cross-referencing right from PseudoJets
1508  if (jetp_index != Invalid) {
1509  assert(jetp_index >= 0);
1510  //cout << _jets.size() <<" "<<jetp_index<<"\n";
1511  _jets[jetp_index].set_cluster_hist_index(local_step);
1512  _set_structure_shared_ptr(_jets[jetp_index]);
1513  }
1514 
1515  if (_writeout_combinations) {
1516  cout << local_step << ": "
1517  << parent1 << " with " << parent2
1518  << "; y = "<< dij<<endl;
1519  }
1520 
1521 }
1522 
1523 
1524 
1525 
1526 //======================================================================
1527 // Return an order in which to read the history such that _history[order[i]]
1528 // will always correspond to the same set of consituent particles if
1529 // two branching histories are equivalent in terms of the particles
1530 // contained in any given pseudojet.
1531 vector<int> ClusterSequence::unique_history_order() const {
1532 
1533  // first construct an array that will tell us the lowest constituent
1534  // of a given jet -- this will always be one of the original
1535  // particles, whose order is well defined and so will help us to
1536  // follow the tree in a unique manner.
1537  valarray<int> lowest_constituent(_history.size());
1538  int hist_n = _history.size();
1539  lowest_constituent = hist_n; // give it a large number
1540  for (int i = 0; i < hist_n; i++) {
1541  // sets things up for the initial partons
1542  lowest_constituent[i] = min(lowest_constituent[i],i);
1543  // propagates them through to the children of this parton
1544  if (_history[i].child > 0) lowest_constituent[_history[i].child]
1545  = min(lowest_constituent[_history[i].child],lowest_constituent[i]);
1546  }
1547 
1548  // establish an array for what we have and have not extracted so far
1549  valarray<bool> extracted(_history.size()); extracted = false;
1550  vector<int> unique_tree;
1551  unique_tree.reserve(_history.size());
1552 
1553  // now work our way through the tree
1554  for (unsigned i = 0; i < n_particles(); i++) {
1555  if (!extracted[i]) {
1556  unique_tree.push_back(i);
1557  extracted[i] = true;
1558  _extract_tree_children(i, extracted, lowest_constituent, unique_tree);
1559  }
1560  }
1561 
1562  return unique_tree;
1563 }
1564 
1565 //======================================================================
1566 // helper for unique_history_order
1567 void ClusterSequence::_extract_tree_children(
1568  int position,
1569  valarray<bool> & extracted,
1570  const valarray<int> & lowest_constituent,
1571  vector<int> & unique_tree) const {
1572  if (!extracted[position]) {
1573  // that means we may have unidentified parents around, so go and
1574  // collect them (extracted[position]) will then be made true)
1575  _extract_tree_parents(position,extracted,lowest_constituent,unique_tree);
1576  }
1577 
1578  // now look after the children...
1579  int child = _history[position].child;
1580  if (child >= 0) _extract_tree_children(child,extracted,lowest_constituent,unique_tree);
1581 }
1582 
1583 
1584 //======================================================================
1585 // return the list of unclustered particles
1586 vector<PseudoJet> ClusterSequence::unclustered_particles() const {
1587  vector<PseudoJet> unclustered;
1588  for (unsigned i = 0; i < n_particles() ; i++) {
1589  if (_history[i].child == Invalid)
1590  unclustered.push_back(_jets[_history[i].jetp_index]);
1591  }
1592  return unclustered;
1593 }
1594 
1595 //======================================================================
1596 /// Return the list of pseudojets in the ClusterSequence that do not
1597 /// have children (and are not among the inclusive jets). They may
1598 /// result from a clustering step or may be one of the pseudojets
1599 /// returned by unclustered_particles().
1600 vector<PseudoJet> ClusterSequence::childless_pseudojets() const {
1601  vector<PseudoJet> unclustered;
1602  for (unsigned i = 0; i < _history.size() ; i++) {
1603  if ((_history[i].child == Invalid) && (_history[i].parent2 != BeamJet))
1604  unclustered.push_back(_jets[_history[i].jetp_index]);
1605  }
1606  return unclustered;
1607 }
1608 
1609 
1610 
1611 //----------------------------------------------------------------------
1612 // returns true if the cluster sequence contains this jet (i.e. jet's
1613 // structure is this cluster sequence's and the cluster history index
1614 // is in a consistent range)
1615 bool ClusterSequence::contains(const PseudoJet & jet) const {
1616  return jet.cluster_hist_index() >= 0
1617  && jet.cluster_hist_index() < int(_history.size())
1619  && jet.associated_cluster_sequence() == this;
1620 }
1621 
1622 
1623 
1624 //======================================================================
1625 // helper for unique_history_order
1626 void ClusterSequence::_extract_tree_parents(
1627  int position,
1628  valarray<bool> & extracted,
1629  const valarray<int> & lowest_constituent,
1630  vector<int> & unique_tree) const {
1631 
1632  if (!extracted[position]) {
1633  int parent1 = _history[position].parent1;
1634  int parent2 = _history[position].parent2;
1635  // where relevant order parents so that we will first treat the
1636  // one containing the smaller "lowest_constituent"
1637  if (parent1 >= 0 && parent2 >= 0) {
1638  if (lowest_constituent[parent1] > lowest_constituent[parent2])
1639  std::swap(parent1, parent2);
1640  }
1641  // then actually run through the parents to extract the constituents...
1642  if (parent1 >= 0 && !extracted[parent1])
1643  _extract_tree_parents(parent1,extracted,lowest_constituent,unique_tree);
1644  if (parent2 >= 0 && !extracted[parent2])
1645  _extract_tree_parents(parent2,extracted,lowest_constituent,unique_tree);
1646  // finally declare this position to be accounted for and push it
1647  // onto our list.
1648  unique_tree.push_back(position);
1649  extracted[position] = true;
1650  }
1651 }
1652 
1653 
1654 //======================================================================
1655 /// carries out the bookkeeping associated with the step of recombining
1656 /// jet_i and jet_j (assuming a distance dij) and returns the index
1657 /// of the recombined jet, newjet_k.
1658 void ClusterSequence::_do_ij_recombination_step(
1659  const int jet_i, const int jet_j,
1660  const double dij,
1661  int & newjet_k) {
1662 
1663  // Create the new jet by recombining the first two.
1664  //
1665  // For efficiency reasons, use a ctr that initialises only the
1666  // shared pointers, since the rest of the info will anyway be dealt
1667  // with by the recombiner.
1668  PseudoJet newjet(false);
1669  _jet_def.recombiner()->recombine(_jets[jet_i], _jets[jet_j], newjet);
1670  _jets.push_back(newjet);
1671  // original version...
1672  //_jets.push_back(_jets[jet_i] + _jets[jet_j]);
1673 
1674  // get its index
1675  newjet_k = _jets.size()-1;
1676 
1677  // get history index
1678  int newstep_k = _history.size();
1679  // and provide jet with the info
1680  _jets[newjet_k].set_cluster_hist_index(newstep_k);
1681 
1682  // finally sort out the history
1683  int hist_i = _jets[jet_i].cluster_hist_index();
1684  int hist_j = _jets[jet_j].cluster_hist_index();
1685 
1686  _add_step_to_history(min(hist_i, hist_j), max(hist_i,hist_j),
1687  newjet_k, dij);
1688 
1689  // _add_step_to_history(newstep_k, min(hist_i, hist_j), max(hist_i,hist_j),
1690  // newjet_k, dij);
1691 
1692 
1693 }
1694 
1695 
1696 //======================================================================
1697 /// carries out the bookkeeping associated with the step of recombining
1698 /// jet_i with the beam
1699 void ClusterSequence::_do_iB_recombination_step(
1700  const int jet_i, const double diB) {
1701  // recombine the jet with the beam
1702  _add_step_to_history(_jets[jet_i].cluster_hist_index(),BeamJet,
1703  Invalid, diB);
1704 
1705  // // get history index
1706  // int newstep_k = _history.size();
1707  //
1708  // _add_step_to_history(newstep_k,_jets[jet_i].cluster_hist_index(),BeamJet,
1709  // Invalid, diB);
1710 
1711 }
1712 
1713 
1714 // make sure the static member _changed_strategy_warning is defined.
1715 LimitedWarning ClusterSequence::_changed_strategy_warning;
1716 
1717 
1718 //----------------------------------------------------------------------
1719 void ClusterSequence::_set_structure_shared_ptr(PseudoJet & j) {
1720  j.set_structure_shared_ptr(_structure_shared_ptr);
1721  // record the use count of the structure shared point to help
1722  // in case we want to ask the CS to handle its own memory
1723  _update_structure_use_count();
1724 }
1725 
1726 
1727 //----------------------------------------------------------------------
1728 void ClusterSequence::_update_structure_use_count() {
1729  // record the use count of the structure shared point to help
1730  // in case we want to ask the CS to handle its own memory
1731  _structure_use_count_after_construction = _structure_shared_ptr.use_count();
1732 }
1733 
1734 //----------------------------------------------------------------------
1735 /// by calling this routine you tell the ClusterSequence to delete
1736 /// itself when all the Pseudojets associated with it have gone out
1737 /// of scope.
1738 void ClusterSequence::delete_self_when_unused() {
1739  // the trick we use to handle this is to modify the use count;
1740  // that way the structure will be deleted when there are no external
1741  // objects left associated the CS and the structure's destructor will then
1742  // look after deleting the cluster sequence
1743 
1744  // first make sure that there is at least one other object
1745  // associated with the CS
1746  int new_count = _structure_shared_ptr.use_count() - _structure_use_count_after_construction;
1747  if (new_count <= 0) {
1748  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");
1749  }
1750 
1751  _structure_shared_ptr.set_count(new_count);
1752  _deletes_self_when_unused = true;
1753 }
1754 
1755 
1756 FASTJET_END_NAMESPACE
1757 
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&#39;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.
the value for the jet algorithm in a JetDefinition for which no algorithm has yet been defined ...
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:772
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&#39;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:774
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
the automatic strategy choice that was being made in FJ 3.0 (restricted to strategies that were prese...
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...
class corresponding to critical internal errors
Definition: Error.hh:109
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&#39;s nearest neighbour (NN) if that particle&#39;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&#39;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