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