Fork me on GitHub

source: svn/trunk/external/fastjet/ClusterSequence.cc@ 1383

Last change on this file since 1383 was 1332, checked in by Pavel Demin, 11 years ago

upgrade fastjet to 3.0.6

  • Property svn:eol-style set to native
  • Property svn:keywords set to Id Revision Date
File size: 56.1 KB
Line 
1//STARTHEADER
2// $Id: ClusterSequence.cc 1332 2013-11-20 20:52:59Z pavel $
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
43FASTJET_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
127using 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.
139std::ostream * ClusterSequence::_fastjet_banner_ostr = &cout;
140
141
142// destructor that guarantees proper bookkeeping for the CS Structure
143ClusterSequence::~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//-----------
169void 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//----------------------------------------------------------------------
199void 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//----------------------------------------------------------------------
211void 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.
360bool ClusterSequence::_first_time = true;
361int ClusterSequence::_n_exclusive_warnings = 0;
362
363
364//----------------------------------------------------------------------
365// the version string
366string fastjet_version_string() {
367 return "FastJet version "+string(fastjet_version);
368}
369
370
371//----------------------------------------------------------------------
372// prints a banner on the first call
373void 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
406void 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
419void 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
438void 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
477string 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
511double 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
587void 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
648void 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
664vector<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
716int 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
735vector<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.
744vector<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
761vector<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
831double 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.
843double 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.
854std::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()
877int 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.
890std::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).
907std::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.
936double 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.
957double 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
981void 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
1015bool 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
1041bool 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
1067bool 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
1089bool 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
1110bool 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
1137vector<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]
1152void 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
1171void 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;
1197vector<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
1228void 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
1258void 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.
1303vector<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
1339void 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
1358vector<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().
1372vector<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)
1387bool ClusterSequence::contains(const PseudoJet & jet) const {
1388 return jet.cluster_hist_index() >= 0
1389 && jet.cluster_hist_index() < int(_history.size())
1390 && jet.has_valid_cluster_sequence()
1391 && jet.associated_cluster_sequence() == this;
1392}
1393
1394
1395
1396//======================================================================
1397// helper for unique_history_order
1398void 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.
1430void 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
1467void 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.
1480LimitedWarning ClusterSequence::_changed_strategy_warning;
1481
1482
1483//----------------------------------------------------------------------
1484void 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//----------------------------------------------------------------------
1493void 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.
1503void 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
1521FASTJET_END_NAMESPACE
1522
Note: See TracBrowser for help on using the repository browser.