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