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 
fastjet::N2Tiled
@ N2Tiled
fastest from about 50..500
Definition: JetDefinition.hh:52
fastjet::NlnNCam2pi2R
@ NlnNCam2pi2R
Chan's closest pair method (in a variant with 2pi+2R coverage), for use exclusively with the Cambridg...
Definition: JetDefinition.hh:76
fastjet::N2MinHeapTiled
@ N2MinHeapTiled
fastest form about 500..10^4
Definition: JetDefinition.hh:50
fastjet::JetDefinition
Definition: JetDefinition.hh:173
fastjet::ClusterSequence::_jets
std::vector< PseudoJet > _jets
This contains the physical PseudoJets; for each PseudoJet one can find the corresponding position in ...
Definition: ClusterSequence.hh:646
fastjet::ClusterSequence::_history
std::vector< history_element > _history
this vector will contain the branching history; for each stage, _history[i].jetp_index indicates wher...
Definition: ClusterSequence.hh:652
fastjet::ClusterSequence
Definition: ClusterSequence.hh:59
fastjet::N2PoorTiled
@ N2PoorTiled
legacy
Definition: JetDefinition.hh:54
fastjet::ClusterSequence::history_element::parent2
int parent2
index in _history where first parent of this jet was created (InexistentParent if this jet is an orig...
Definition: ClusterSequence.hh:409
fastjet::ClusterSequence::history_element::jetp_index
int jetp_index
index in _history where the current jet is recombined with another jet to form its child.
Definition: ClusterSequence.hh:420
fastjet::cambridge_algorithm
@ cambridge_algorithm
the longitudinally invariant variant of the cambridge algorithm (aka Aachen algoithm).
Definition: JetDefinition.hh:94
fastjet::NlnNCam4pi
@ NlnNCam4pi
Chan's closest pair method (in a variant with 4pi coverage), for use exclusively with the Cambridge a...
Definition: JetDefinition.hh:72
fastjet::PseudoJet::has_valid_cluster_sequence
bool has_valid_cluster_sequence() const
returns true if this PseudoJet has an associated and still valid(ated) ClusterSequence.
Definition: PseudoJet.cc:433
fastjet::ClusterSequenceStructure::set_associated_cs
virtual void set_associated_cs(const ClusterSequence *new_cs)
set the associated csw
Definition: ClusterSequenceStructure.hh:105
fastjet::PseudoJet::associated_cluster_sequence
const ClusterSequence * associated_cluster_sequence() const
get a (const) pointer to the parent ClusterSequence (NULL if inexistent)
Definition: PseudoJet.cc:423
fastjet::fastjet_version_string
string fastjet_version_string()
return a string containing information about the release
Definition: ClusterSequence.cc:366
fastjet::ClusterSequence::history_element
Definition: ClusterSequence.hh:404
fastjet::FunctionOfPseudoJet< PseudoJet >
fastjet::NlnN
@ NlnN
best of the NlnN variants – best overall for N>10^4.
Definition: JetDefinition.hh:63
fastjet::PseudoJet::set_structure_shared_ptr
void set_structure_shared_ptr(const SharedPtr< PseudoJetStructureBase > &structure)
set the associated structure
Definition: PseudoJet.cc:450
fastjet::genkt_algorithm
@ genkt_algorithm
like the k_t but with distance measures dij = min(kti^{2p},ktj^{2p}) Delta R_{ij}^2 / R^2 diB = 1/kti...
Definition: JetDefinition.hh:103
fastjet::PseudoJet::kt2
double kt2() const
returns the squared transverse momentum
Definition: PseudoJet.hh:143
fastjet::PseudoJet::cluster_hist_index
int cluster_hist_index() const
return the cluster_hist_index, intended to be used by clustering routines.
Definition: PseudoJet.hh:764
fastjet::antikt_algorithm
@ antikt_algorithm
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
Definition: JetDefinition.hh:98
fastjet::PseudoJet
Definition: PseudoJet.hh:65
fastjet::ClusterSequence::history_element::dij
double dij
index in the _jets vector where we will find the
Definition: ClusterSequence.hh:427
fastjet::ee_genkt_algorithm
@ ee_genkt_algorithm
the e+e- genkt algorithm (R > 2 and p=1 gives ee_kt)
Definition: JetDefinition.hh:114
fastjet::NlnN4pi
@ NlnN4pi
legacy N ln N using 4pi coverage of cylinder
Definition: JetDefinition.hh:68
fastjet::NlnN3pi
@ NlnN3pi
legacy N ln N using 3pi coverage of cylinder.
Definition: JetDefinition.hh:66
fastjet::PseudoJet::set_cluster_hist_index
void set_cluster_hist_index(const int index)
set the cluster_hist_index, intended to be used by clustering routines.
Definition: PseudoJet.hh:766
fastjet::Strategy
Strategy
the various options for the algorithmic strategy to adopt in clustering events with kt and cambridge ...
Definition: JetDefinition.hh:48
fastjet::LimitedWarning
Definition: LimitedWarning.hh:45
fastjet::ClusterSequence::history_element::max_dij_so_far
double max_dij_so_far
the distance corresponding to the recombination at this stage of the clustering.
Definition: ClusterSequence.hh:430
fastjet::ClusterSequence::history_element::child
int child
index in _history where second parent of this jet was created (InexistentParent if this jet is an ori...
Definition: ClusterSequence.hh:415
fastjet::NlnNCam
@ NlnNCam
Chan's closest pair method (in a variant with 2pi+minimal extra variant), for use exclusively with th...
Definition: JetDefinition.hh:80
fastjet::plugin_strategy
@ plugin_strategy
the plugin has been used...
Definition: JetDefinition.hh:82
fastjet::N2Plain
@ N2Plain
fastest below 50
Definition: JetDefinition.hh:56
fastjet::PseudoJet::perp2
double perp2() const
returns the squared transverse momentum
Definition: PseudoJet.hh:139
fastjet::N3Dumb
@ N3Dumb
worse even than the usual N^3 algorithms
Definition: JetDefinition.hh:58
fastjet::plugin_algorithm
@ plugin_algorithm
any plugin algorithm supplied by the user
Definition: JetDefinition.hh:117
fastjet::ClusterSequenceStructure
Definition: ClusterSequenceStructure.hh:59
fastjet::ee_kt_algorithm
@ ee_kt_algorithm
the e+e- kt algorithm
Definition: JetDefinition.hh:112
fastjet::Error
Definition: Error.hh:41
fastjet::cambridge_for_passive_algorithm
@ cambridge_for_passive_algorithm
a version of cambridge with a special distance measure for particles whose pt is < extra_param()
Definition: JetDefinition.hh:106
fastjet::kt_algorithm
@ kt_algorithm
the longitudinally invariant kt algorithm
Definition: JetDefinition.hh:91