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