Fork me on GitHub

source: git/external/fastjet/ClusterSequence.cc@ 49234af

ImprovedOutputFile Timing dual_readout llp
Last change on this file since 49234af was 273e668, checked in by Pavel Demin <pavel.demin@…>, 10 years ago

upgrade FastJet to version 3.1.0

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