Fork me on GitHub

source: git/external/fastjet/ClusterSequence.cc@ 947209d

ImprovedOutputFile Timing dual_readout llp
Last change on this file since 947209d was 1d208a2, checked in by Pavel Demin <pavel.demin@…>, 8 years ago

update FastJet library to 3.2.1 and Nsubjettiness library to 2.2.4

  • Property mode set to 100644
File size: 66.7 KB
Line 
1//FJSTARTHEADER
2// $Id: ClusterSequence.cc 4154 2016-07-20 16:20:48Z soyez $
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.get());
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 } else if (_strategy == N2MHTLazy9AntiKtSeparateGhosts) {
375#ifndef __FJCORE__
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
729 assert(0 && "Code should never reach here");
730
731 return N2MHTLazy9;
732
733}
734
735
736// //----------------------------------------------------------------------
737// /// transfer the sequence contained in other_seq into our own;
738// /// any plugin "extras" contained in the from_seq will be lost
739// /// from there.
740// void ClusterSequence::transfer_from_sequence(ClusterSequence & from_seq) {
741//
742// if (will_delete_self_when_unused())
743// throw(Error("cannot use CS::transfer_from_sequence after a call to delete_self_when_unused()"));
744//
745// // the metadata
746// _jet_def = from_seq._jet_def ;
747// _writeout_combinations = from_seq._writeout_combinations ;
748// _initial_n = from_seq._initial_n ;
749// _Rparam = from_seq._Rparam ;
750// _R2 = from_seq._R2 ;
751// _invR2 = from_seq._invR2 ;
752// _strategy = from_seq._strategy ;
753// _jet_algorithm = from_seq._jet_algorithm ;
754// _plugin_activated = from_seq._plugin_activated ;
755//
756// // the data
757// _jets = from_seq._jets;
758// _history = from_seq._history;
759// // the following transfers ownership of the extras from the from_seq
760// _extras = from_seq._extras;
761//
762// // transfer of ownership
763// if (_structure_shared_ptr()) {
764// // anything that is currently associated with the cluster sequence
765// // should be told that its cluster sequence no longer exists
766// ClusterSequenceStructure* csi = dynamic_cast<ClusterSequenceStructure*>(_structure_shared_ptr());
767// assert(csi != NULL);
768// csi->set_associated_cs(NULL);
769// }
770// // create a new _structure_shared_ptr to reflect the fact that
771// // this CS is essentially a new one
772// _structure_shared_ptr.reset(new ClusterSequenceStructure(this));
773// _update_structure_use_count();
774//
775// for (vector<PseudoJet>::iterator jit = _jets.begin(); jit != _jets.end(); jit++)
776// _set_structure_shared_ptr(*jit);
777// }
778
779
780ClusterSequence & ClusterSequence::operator=(const ClusterSequence & cs) {
781 // self assignment is trivial
782 if (&cs != this) {
783 _deletes_self_when_unused = false;
784 transfer_from_sequence(cs);
785 }
786 return *this;
787}
788
789//----------------------------------------------------------------------
790// transfer the sequence contained in other_seq into our own;
791// any plugin "extras" contained in the from_seq will be lost
792// from there.
793//
794// It also sets the ClusterSequence pointers of the PseudoJets in
795// the history to point to this ClusterSequence
796//
797// The second argument is an action that will be applied on every
798// jets in the resulting ClusterSequence
799void ClusterSequence::transfer_from_sequence(const ClusterSequence & from_seq,
800 const FunctionOfPseudoJet<PseudoJet> * action_on_jets){
801
802 if (will_delete_self_when_unused())
803 throw(Error("cannot use CS::transfer_from_sequence after a call to delete_self_when_unused()"));
804
805 // the metadata
806 _jet_def = from_seq._jet_def ;
807 _writeout_combinations = from_seq._writeout_combinations ;
808 _initial_n = from_seq._initial_n ;
809 _Rparam = from_seq._Rparam ;
810 _R2 = from_seq._R2 ;
811 _invR2 = from_seq._invR2 ;
812 _strategy = from_seq._strategy ;
813 _jet_algorithm = from_seq._jet_algorithm ;
814 _plugin_activated = from_seq._plugin_activated ;
815
816 // the data
817
818 // apply the transformation on the jets if needed
819 if (action_on_jets)
820 _jets = (*action_on_jets)(from_seq._jets);
821 else
822 _jets = from_seq._jets;
823 _history = from_seq._history;
824 // the following shares ownership of the extras with the from_seq;
825 // no transformations will be applied to the extras
826 _extras = from_seq._extras;
827
828 // clean up existing structure
829 if (_structure_shared_ptr) {
830 // If there are jets associated with an old version of the CS and
831 // a new one, keeping track of when to delete the CS becomes more
832 // complex; so we don't allow this situation to occur.
833 if (_deletes_self_when_unused) throw Error("transfer_from_sequence cannot be used for a cluster sequence that deletes self when unused");
834
835 // anything that is currently associated with the cluster sequence
836 // should be told that its cluster sequence no longer exists
837 ClusterSequenceStructure* csi = dynamic_cast<ClusterSequenceStructure*>(_structure_shared_ptr.get());
838 assert(csi != NULL);
839 csi->set_associated_cs(NULL);
840 }
841 // create a new _structure_shared_ptr to reflect the fact that
842 // this CS is essentially a new one
843 _structure_shared_ptr.reset(new ClusterSequenceStructure(this));
844 _update_structure_use_count();
845
846 for (unsigned int i=0; i<_jets.size(); i++){
847 // we reset the cluster history index in case action_on_jets
848 // messed up with it
849 _jets[i].set_cluster_hist_index(from_seq._jets[i].cluster_hist_index());
850
851 // reset the structure pointer
852 _set_structure_shared_ptr(_jets[i]);
853 }
854}
855
856
857//----------------------------------------------------------------------
858// record an ij recombination and reset the _jets[newjet_k] momentum and
859// user index to be those of newjet
860void ClusterSequence::plugin_record_ij_recombination(
861 int jet_i, int jet_j, double dij,
862 const PseudoJet & newjet, int & newjet_k) {
863
864 plugin_record_ij_recombination(jet_i, jet_j, dij, newjet_k);
865
866 // now transfer newjet into place
867 int tmp_index = _jets[newjet_k].cluster_hist_index();
868 _jets[newjet_k] = newjet;
869 _jets[newjet_k].set_cluster_hist_index(tmp_index);
870 _set_structure_shared_ptr(_jets[newjet_k]);
871}
872
873
874//----------------------------------------------------------------------
875// return all inclusive jets with pt > ptmin
876vector<PseudoJet> ClusterSequence::inclusive_jets (const double ptmin) const{
877 double dcut = ptmin*ptmin;
878 int i = _history.size() - 1; // last jet
879 vector<PseudoJet> jets_local;
880 if (_jet_algorithm == kt_algorithm) {
881 while (i >= 0) {
882 // with our specific definition of dij and diB (i.e. R appears only in
883 // dij), then dij==diB is the same as the jet.perp2() and we can exploit
884 // this in selecting the jets...
885 if (_history[i].max_dij_so_far < dcut) {break;}
886 if (_history[i].parent2 == BeamJet && _history[i].dij >= dcut) {
887 // for beam jets
888 int parent1 = _history[i].parent1;
889 jets_local.push_back(_jets[_history[parent1].jetp_index]);}
890 i--;
891 }
892 } else if (_jet_algorithm == cambridge_algorithm) {
893 while (i >= 0) {
894 // inclusive jets are all at end of clustering sequence in the
895 // Cambridge algorithm -- so if we find a non-exclusive jet, then
896 // we can exit
897 if (_history[i].parent2 != BeamJet) {break;}
898 int parent1 = _history[i].parent1;
899 const PseudoJet & jet = _jets[_history[parent1].jetp_index];
900 if (jet.perp2() >= dcut) {jets_local.push_back(jet);}
901 i--;
902 }
903 } else if (_jet_algorithm == plugin_algorithm
904 || _jet_algorithm == ee_kt_algorithm
905 || _jet_algorithm == antikt_algorithm
906 || _jet_algorithm == genkt_algorithm
907 || _jet_algorithm == ee_genkt_algorithm
908 || _jet_algorithm == cambridge_for_passive_algorithm) {
909 // for inclusive jets with a plugin algorithm, we make no
910 // assumptions about anything (relation of dij to momenta,
911 // ordering of the dij, etc.)
912 while (i >= 0) {
913 if (_history[i].parent2 == BeamJet) {
914 int parent1 = _history[i].parent1;
915 const PseudoJet & jet = _jets[_history[parent1].jetp_index];
916 if (jet.perp2() >= dcut) {jets_local.push_back(jet);}
917 }
918 i--;
919 }
920 } else {throw Error("cs::inclusive_jets(...): Unrecognized jet algorithm");}
921 return jets_local;
922}
923
924
925//----------------------------------------------------------------------
926// return the number of exclusive jets that would have been obtained
927// running the algorithm in exclusive mode with the given dcut
928int ClusterSequence::n_exclusive_jets (const double dcut) const {
929
930 // first locate the point where clustering would have stopped (i.e. the
931 // first time max_dij_so_far > dcut)
932 int i = _history.size() - 1; // last jet
933 while (i >= 0) {
934 if (_history[i].max_dij_so_far <= dcut) {break;}
935 i--;
936 }
937 int stop_point = i + 1;
938 // relation between stop_point, njets assumes one extra jet disappears
939 // at each clustering.
940 int njets = 2*_initial_n - stop_point;
941 return njets;
942}
943
944//----------------------------------------------------------------------
945// return all exclusive jets that would have been obtained running
946// the algorithm in exclusive mode with the given dcut
947vector<PseudoJet> ClusterSequence::exclusive_jets (const double dcut) const {
948 int njets = n_exclusive_jets(dcut);
949 return exclusive_jets(njets);
950}
951
952
953//----------------------------------------------------------------------
954// return the jets obtained by clustering the event to n jets.
955// Throw an error if there are fewer than n particles.
956vector<PseudoJet> ClusterSequence::exclusive_jets (const int njets) const {
957
958 // make sure the user does not ask for more than jets than there
959 // were particles in the first place.
960 if (njets > _initial_n) {
961 ostringstream err;
962 err << "Requested " << njets << " exclusive jets, but there were only "
963 << _initial_n << " particles in the event";
964 throw Error(err.str());
965 }
966
967 return exclusive_jets_up_to(njets);
968}
969
970//----------------------------------------------------------------------
971// return the jets obtained by clustering the event to n jets.
972// If there are fewer than n particles, simply return all particles
973vector<PseudoJet> ClusterSequence::exclusive_jets_up_to (const int njets) const {
974
975 // provide a warning when extracting exclusive jets for algorithms
976 // that does not support it explicitly.
977 // Native algorithm that support it are: kt, ee_kt, Cambridge/Aachen,
978 // genkt and ee_genkt (both with p>=0)
979 // For plugins, we check Plugin::exclusive_sequence_meaningful()
980 if (( _jet_def.jet_algorithm() != kt_algorithm) &&
981 ( _jet_def.jet_algorithm() != cambridge_algorithm) &&
982 ( _jet_def.jet_algorithm() != ee_kt_algorithm) &&
983 (((_jet_def.jet_algorithm() != genkt_algorithm) &&
984 (_jet_def.jet_algorithm() != ee_genkt_algorithm)) ||
985 (_jet_def.extra_param() <0)) &&
986 ((_jet_def.jet_algorithm() != plugin_algorithm) ||
987 (!_jet_def.plugin()->exclusive_sequence_meaningful()))) {
988 _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.");
989 }
990
991
992 // calculate the point where we have to stop the clustering.
993 // relation between stop_point, njets assumes one extra jet disappears
994 // at each clustering.
995 int stop_point = 2*_initial_n - njets;
996 // make sure it's safe when more jets are requested than there are particles
997 if (stop_point < _initial_n) stop_point = _initial_n;
998
999 // some sanity checking to make sure that e+e- does not give us
1000 // surprises (should we ever implement e+e-)...
1001 if (2*_initial_n != static_cast<int>(_history.size())) {
1002 ostringstream err;
1003 err << "2*_initial_n != _history.size() -- this endangers internal assumptions!\n";
1004 throw Error(err.str());
1005 //assert(false);
1006 }
1007
1008 // now go forwards and reconstitute the jets that we have --
1009 // basically for any history element, see if the parent jets to
1010 // which it refers were created before the stopping point -- if they
1011 // were then add them to the list, otherwise they are subsequent
1012 // recombinations of the jets that we are looking for.
1013 vector<PseudoJet> jets_local;
1014 for (unsigned int i = stop_point; i < _history.size(); i++) {
1015 int parent1 = _history[i].parent1;
1016 if (parent1 < stop_point) {
1017 jets_local.push_back(_jets[_history[parent1].jetp_index]);
1018 }
1019 int parent2 = _history[i].parent2;
1020 if (parent2 < stop_point && parent2 > 0) {
1021 jets_local.push_back(_jets[_history[parent2].jetp_index]);
1022 }
1023
1024 }
1025
1026 // sanity check...
1027 if (int(jets_local.size()) != min(_initial_n, njets)) {
1028 ostringstream err;
1029 err << "ClusterSequence::exclusive_jets: size of returned vector ("
1030 <<jets_local.size()<<") does not coincide with requested number of jets ("
1031 <<njets<<")";
1032 throw Error(err.str());
1033 }
1034
1035 return jets_local;
1036}
1037
1038//----------------------------------------------------------------------
1039/// return the dmin corresponding to the recombination that went from
1040/// n+1 to n jets
1041double ClusterSequence::exclusive_dmerge (const int njets) const {
1042 assert(njets >= 0);
1043 if (njets >= _initial_n) {return 0.0;}
1044 return _history[2*_initial_n-njets-1].dij;
1045}
1046
1047
1048//----------------------------------------------------------------------
1049/// return the maximum of the dmin encountered during all recombinations
1050/// up to the one that led to an n-jet final state; identical to
1051/// exclusive_dmerge, except in cases where the dmin do not increase
1052/// monotonically.
1053double ClusterSequence::exclusive_dmerge_max (const int njets) const {
1054 assert(njets >= 0);
1055 if (njets >= _initial_n) {return 0.0;}
1056 return _history[2*_initial_n-njets-1].max_dij_so_far;
1057}
1058
1059
1060//----------------------------------------------------------------------
1061/// return a vector of all subjets of the current jet (in the sense
1062/// of the exclusive algorithm) that would be obtained when running
1063/// the algorithm with the given dcut.
1064std::vector<PseudoJet> ClusterSequence::exclusive_subjets
1065 (const PseudoJet & jet, const double dcut) const {
1066
1067 set<const history_element*> subhist;
1068
1069 // get the set of history elements that correspond to subjets at
1070 // scale dcut
1071 get_subhist_set(subhist, jet, dcut, 0);
1072
1073 // now transfer this into a sequence of jets
1074 vector<PseudoJet> subjets;
1075 subjets.reserve(subhist.size());
1076 for (set<const history_element*>::iterator elem = subhist.begin();
1077 elem != subhist.end(); elem++) {
1078 subjets.push_back(_jets[(*elem)->jetp_index]);
1079 }
1080 return subjets;
1081}
1082
1083//----------------------------------------------------------------------
1084/// return the size of exclusive_subjets(...); still n ln n with same
1085/// coefficient, but marginally more efficient than manually taking
1086/// exclusive_subjets.size()
1087int ClusterSequence::n_exclusive_subjets(const PseudoJet & jet,
1088 const double dcut) const {
1089 set<const history_element*> subhist;
1090 // get the set of history elements that correspond to subjets at
1091 // scale dcut
1092 get_subhist_set(subhist, jet, dcut, 0);
1093 return subhist.size();
1094}
1095
1096//----------------------------------------------------------------------
1097/// return the list of subjets obtained by unclustering the supplied
1098/// jet down to nsub subjets. Throws an error if there are fewer than
1099/// nsub particles in the jet.
1100std::vector<PseudoJet> ClusterSequence::exclusive_subjets
1101 (const PseudoJet & jet, int nsub) const {
1102 vector<PseudoJet> subjets = exclusive_subjets_up_to(jet, nsub);
1103 if (int(subjets.size()) < nsub) {
1104 ostringstream err;
1105 err << "Requested " << nsub << " exclusive subjets, but there were only "
1106 << subjets.size() << " particles in the jet";
1107 throw Error(err.str());
1108 }
1109 return subjets;
1110
1111}
1112
1113//----------------------------------------------------------------------
1114/// return the list of subjets obtained by unclustering the supplied
1115/// jet down to nsub subjets (or all constituents if there are fewer
1116/// than nsub).
1117std::vector<PseudoJet> ClusterSequence::exclusive_subjets_up_to
1118 (const PseudoJet & jet, int nsub) const {
1119
1120 set<const history_element*> subhist;
1121
1122 // prepare the vector into which we'll put the result
1123 vector<PseudoJet> subjets;
1124 if (nsub < 0) throw Error("Requested a negative number of subjets. This is nonsensical.");
1125 if (nsub == 0) return subjets;
1126
1127 // get the set of history elements that correspond to subjets at
1128 // scale dcut
1129 get_subhist_set(subhist, jet, -1.0, nsub);
1130
1131 // now transfer this into a sequence of jets
1132 subjets.reserve(subhist.size());
1133 for (set<const history_element*>::iterator elem = subhist.begin();
1134 elem != subhist.end(); elem++) {
1135 subjets.push_back(_jets[(*elem)->jetp_index]);
1136 }
1137 return subjets;
1138}
1139
1140
1141//----------------------------------------------------------------------
1142/// return the dij that was present in the merging nsub+1 -> nsub
1143/// subjets inside this jet.
1144///
1145/// If the jet has nsub or fewer constituents, it will return 0.
1146double ClusterSequence::exclusive_subdmerge(const PseudoJet & jet, int nsub) const {
1147 set<const history_element*> subhist;
1148
1149 // get the set of history elements that correspond to subjets at
1150 // scale dcut
1151 get_subhist_set(subhist, jet, -1.0, nsub);
1152
1153 set<const history_element*>::iterator highest = subhist.end();
1154 highest--;
1155 /// will be zero if nconst <= nsub, since highest will be an original
1156 /// particle have zero dij
1157 return (*highest)->dij;
1158}
1159
1160
1161//----------------------------------------------------------------------
1162/// return the maximum dij that occurred in the whole event at the
1163/// stage that the nsub+1 -> nsub merge of subjets occurred inside
1164/// this jet.
1165///
1166/// If the jet has nsub or fewer constituents, it will return 0.
1167double ClusterSequence::exclusive_subdmerge_max(const PseudoJet & jet, int nsub) const {
1168
1169 set<const history_element*> subhist;
1170
1171 // get the set of history elements that correspond to subjets at
1172 // scale dcut
1173 get_subhist_set(subhist, jet, -1.0, nsub);
1174
1175 set<const history_element*>::iterator highest = subhist.end();
1176 highest--;
1177 /// will be zero if nconst <= nsub, since highest will be an original
1178 /// particle have zero dij
1179 return (*highest)->max_dij_so_far;
1180}
1181
1182
1183
1184//----------------------------------------------------------------------
1185/// return a set of pointers to history entries corresponding to the
1186/// subjets of this jet; one stops going working down through the
1187/// subjets either when
1188/// - there is no further to go
1189/// - one has found maxjet entries
1190/// - max_dij_so_far <= dcut
1191void ClusterSequence::get_subhist_set(set<const history_element*> & subhist,
1192 const PseudoJet & jet,
1193 double dcut, int maxjet) const {
1194 assert(contains(jet));
1195
1196 subhist.clear();
1197 subhist.insert(&(_history[jet.cluster_hist_index()]));
1198
1199 // establish the set of jets that are relevant
1200 int njet = 1;
1201 while (true) {
1202 // first find out if we need to probe deeper into jet.
1203 // Get history element closest to end of sequence
1204 set<const history_element*>::iterator highest = subhist.end();
1205 assert (highest != subhist.begin());
1206 highest--;
1207 const history_element* elem = *highest;
1208 // make sure we haven't got too many jets
1209 if (njet == maxjet) break;
1210 // make sure it has parents
1211 if (elem->parent1 < 0) break;
1212 // make sure that we still resolve it at scale dcut
1213 if (elem->max_dij_so_far <= dcut) break;
1214
1215 // then do so: replace "highest" with its two parents
1216 subhist.erase(highest);
1217 subhist.insert(&(_history[elem->parent1]));
1218 subhist.insert(&(_history[elem->parent2]));
1219 njet++;
1220 }
1221}
1222
1223//----------------------------------------------------------------------
1224// work through the object's history until
1225bool ClusterSequence::object_in_jet(const PseudoJet & object,
1226 const PseudoJet & jet) const {
1227
1228 // make sure the object conceivably belongs to this clustering
1229 // sequence
1230 assert(contains(object) && contains(jet));
1231
1232 const PseudoJet * this_object = &object;
1233 const PseudoJet * childp;
1234 while(true) {
1235 if (this_object->cluster_hist_index() == jet.cluster_hist_index()) {
1236 return true;
1237 } else if (has_child(*this_object, childp)) {
1238 this_object = childp;
1239 } else {
1240 return false;
1241 }
1242 }
1243}
1244
1245//----------------------------------------------------------------------
1246/// if the jet has parents in the clustering, it returns true
1247/// and sets parent1 and parent2 equal to them.
1248///
1249/// if it has no parents it returns false and sets parent1 and
1250/// parent2 to zero
1251bool ClusterSequence::has_parents(const PseudoJet & jet, PseudoJet & parent1,
1252 PseudoJet & parent2) const {
1253
1254 const history_element & hist = _history[jet.cluster_hist_index()];
1255
1256 // make sure we do not run into any unexpected situations --
1257 // i.e. both parents valid, or neither
1258 assert ((hist.parent1 >= 0 && hist.parent2 >= 0) ||
1259 (hist.parent1 < 0 && hist.parent2 < 0));
1260
1261 if (hist.parent1 < 0) {
1262 parent1 = PseudoJet(0.0,0.0,0.0,0.0);
1263 parent2 = parent1;
1264 return false;
1265 } else {
1266 parent1 = _jets[_history[hist.parent1].jetp_index];
1267 parent2 = _jets[_history[hist.parent2].jetp_index];
1268 // order the parents in decreasing pt
1269 if (parent1.perp2() < parent2.perp2()) std::swap(parent1,parent2);
1270 return true;
1271 }
1272}
1273
1274//----------------------------------------------------------------------
1275/// if the jet has a child then return true and give the child jet
1276/// otherwise return false and set the child to zero
1277bool ClusterSequence::has_child(const PseudoJet & jet, PseudoJet & child) const {
1278
1279 //const history_element & hist = _history[jet.cluster_hist_index()];
1280 //
1281 //if (hist.child >= 0) {
1282 // child = _jets[_history[hist.child].jetp_index];
1283 // return true;
1284 //} else {
1285 // child = PseudoJet(0.0,0.0,0.0,0.0);
1286 // return false;
1287 //}
1288 const PseudoJet * childp;
1289 bool res = has_child(jet, childp);
1290 if (res) {
1291 child = *childp;
1292 return true;
1293 } else {
1294 child = PseudoJet(0.0,0.0,0.0,0.0);
1295 return false;
1296 }
1297}
1298
1299bool ClusterSequence::has_child(const PseudoJet & jet, const PseudoJet * & childp) const {
1300
1301 const history_element & hist = _history[jet.cluster_hist_index()];
1302
1303 // check that this jet has a child and that the child corresponds to
1304 // a true jet [RETHINK-IF-CHANGE-NUMBERING: what is the right
1305 // behaviour if the child is the same jet but made inclusive...?]
1306 if (hist.child >= 0 && _history[hist.child].jetp_index >= 0) {
1307 childp = &(_jets[_history[hist.child].jetp_index]);
1308 return true;
1309 } else {
1310 childp = NULL;
1311 return false;
1312 }
1313}
1314
1315
1316//----------------------------------------------------------------------
1317/// if this jet has a child (and so a partner) return true
1318/// and give the partner, otherwise return false and set the
1319/// partner to zero
1320bool ClusterSequence::has_partner(const PseudoJet & jet,
1321 PseudoJet & partner) const {
1322
1323 const history_element & hist = _history[jet.cluster_hist_index()];
1324
1325 // make sure we have a child and that the child does not correspond
1326 // to a clustering with the beam (or some other invalid quantity)
1327 if (hist.child >= 0 && _history[hist.child].parent2 >= 0) {
1328 const history_element & child_hist = _history[hist.child];
1329 if (child_hist.parent1 == jet.cluster_hist_index()) {
1330 // partner will be child's parent2 -- for iB clustering
1331 // parent2 will not be valid
1332 partner = _jets[_history[child_hist.parent2].jetp_index];
1333 } else {
1334 // partner will be child's parent1
1335 partner = _jets[_history[child_hist.parent1].jetp_index];
1336 }
1337 return true;
1338 } else {
1339 partner = PseudoJet(0.0,0.0,0.0,0.0);
1340 return false;
1341 }
1342}
1343
1344
1345//----------------------------------------------------------------------
1346// return a vector of the particles that make up a jet
1347vector<PseudoJet> ClusterSequence::constituents (const PseudoJet & jet) const {
1348 vector<PseudoJet> subjets;
1349 add_constituents(jet, subjets);
1350 return subjets;
1351}
1352
1353//----------------------------------------------------------------------
1354/// output the supplied vector of jets in a format that can be read
1355/// by an appropriate root script; the format is:
1356/// jet-n jet-px jet-py jet-pz jet-E
1357/// particle-n particle-rap particle-phi particle-pt
1358/// particle-n particle-rap particle-phi particle-pt
1359/// ...
1360/// #END
1361/// ... [i.e. above repeated]
1362void ClusterSequence::print_jets_for_root(const std::vector<PseudoJet> & jets_in,
1363 ostream & ostr) const {
1364 for (unsigned i = 0; i < jets_in.size(); i++) {
1365 ostr << i << " "
1366 << jets_in[i].px() << " "
1367 << jets_in[i].py() << " "
1368 << jets_in[i].pz() << " "
1369 << jets_in[i].E() << endl;
1370 vector<PseudoJet> cst = constituents(jets_in[i]);
1371 for (unsigned j = 0; j < cst.size() ; j++) {
1372 ostr << " " << j << " "
1373 << cst[j].rap() << " "
1374 << cst[j].phi() << " "
1375 << cst[j].perp() << endl;
1376 }
1377 ostr << "#END" << endl;
1378 }
1379}
1380
1381void ClusterSequence::print_jets_for_root(const std::vector<PseudoJet> & jets_in,
1382 const std::string & filename,
1383 const std::string & comment ) const {
1384 std::ofstream ostr(filename.c_str());
1385 if (comment != "") ostr << "# " << comment << endl;
1386 print_jets_for_root(jets_in, ostr);
1387}
1388
1389
1390// Not yet. Perhaps in a future release
1391// //----------------------------------------------------------------------
1392// // print out all inclusive jets with pt > ptmin
1393// void ClusterSequence::print_jets (const double ptmin) const{
1394// vector<PseudoJet> jets = sorted_by_pt(inclusive_jets(ptmin));
1395//
1396// for (size_t j = 0; j < jets.size(); j++) {
1397// printf("%5u %7.3f %7.3f %9.3f\n",
1398// j,jets[j].rap(),jets[j].phi(),jets[j].perp());
1399// }
1400// }
1401
1402//----------------------------------------------------------------------
1403/// returns a vector of size n_particles() which indicates, for
1404/// each of the initial particles (in the order in which they were
1405/// supplied), which of the supplied jets it belongs to; if it does
1406/// not belong to any of the supplied jets, the index is set to -1;
1407vector<int> ClusterSequence::particle_jet_indices(
1408 const vector<PseudoJet> & jets_in) const {
1409
1410 vector<int> indices(n_particles());
1411
1412 // first label all particles as not belonging to any jets
1413 for (unsigned ipart = 0; ipart < n_particles(); ipart++)
1414 indices[ipart] = -1;
1415
1416 // then for each of the jets relabel its consituents as belonging to
1417 // that jet
1418 for (unsigned ijet = 0; ijet < jets_in.size(); ijet++) {
1419
1420 vector<PseudoJet> jet_constituents(constituents(jets_in[ijet]));
1421
1422 for (unsigned ip = 0; ip < jet_constituents.size(); ip++) {
1423 // a safe (if slightly redundant) way of getting the particle
1424 // index (for initial particles it is actually safe to assume
1425 // ipart=iclust).
1426 unsigned iclust = jet_constituents[ip].cluster_hist_index();
1427 unsigned ipart = history()[iclust].jetp_index;
1428 indices[ipart] = ijet;
1429 }
1430 }
1431
1432 return indices;
1433}
1434
1435
1436//----------------------------------------------------------------------
1437// recursive routine that adds on constituents of jet to the subjet_vector
1438void ClusterSequence::add_constituents (
1439 const PseudoJet & jet, vector<PseudoJet> & subjet_vector) const {
1440 // find out position in cluster history
1441 int i = jet.cluster_hist_index();
1442 int parent1 = _history[i].parent1;
1443 int parent2 = _history[i].parent2;
1444
1445 if (parent1 == InexistentParent) {
1446 // It is an original particle (labelled by its parent having value
1447 // InexistentParent), therefore add it on to the subjet vector
1448 // Note: we add the initial particle and not simply 'jet' so that
1449 // calling add_constituents with a subtracted jet containing
1450 // only one particle will work.
1451 subjet_vector.push_back(_jets[i]);
1452 return;
1453 }
1454
1455 // add parent 1
1456 add_constituents(_jets[_history[parent1].jetp_index], subjet_vector);
1457
1458 // see if parent2 is a real jet; if it is then add its constituents
1459 if (parent2 != BeamJet) {
1460 add_constituents(_jets[_history[parent2].jetp_index], subjet_vector);
1461 }
1462}
1463
1464
1465
1466//----------------------------------------------------------------------
1467// initialise the history in a standard way
1468void ClusterSequence::_add_step_to_history (
1469 //NO_LONGER_USED: const int step_number,
1470 const int parent1,
1471 const int parent2, const int jetp_index,
1472 const double dij) {
1473
1474 history_element element;
1475 element.parent1 = parent1;
1476 element.parent2 = parent2;
1477 element.jetp_index = jetp_index;
1478 element.child = Invalid;
1479 element.dij = dij;
1480 element.max_dij_so_far = max(dij,_history[_history.size()-1].max_dij_so_far);
1481 _history.push_back(element);
1482
1483 int local_step = _history.size()-1;
1484 //#ifndef __NO_ASSERTS__
1485 //assert(local_step == step_number);
1486 //#endif
1487
1488 // sanity check: make sure the particles have not already been recombined
1489 //
1490 // Note that good practice would make this an assert (since this is
1491 // a serious internal issue). However, we decided to throw an
1492 // InternalError so that the end user can decide to catch it and
1493 // retry the clustering with a different strategy.
1494
1495 assert(parent1 >= 0);
1496 if (_history[parent1].child != Invalid){
1497 throw InternalError("trying to recomine an object that has previsously been recombined");
1498 }
1499 _history[parent1].child = local_step;
1500 if (parent2 >= 0) {
1501 if (_history[parent2].child != Invalid){
1502 throw InternalError("trying to recomine an object that has previsously been recombined");
1503 }
1504 _history[parent2].child = local_step;
1505 }
1506
1507 // get cross-referencing right from PseudoJets
1508 if (jetp_index != Invalid) {
1509 assert(jetp_index >= 0);
1510 //cout << _jets.size() <<" "<<jetp_index<<"\n";
1511 _jets[jetp_index].set_cluster_hist_index(local_step);
1512 _set_structure_shared_ptr(_jets[jetp_index]);
1513 }
1514
1515 if (_writeout_combinations) {
1516 cout << local_step << ": "
1517 << parent1 << " with " << parent2
1518 << "; y = "<< dij<<endl;
1519 }
1520
1521}
1522
1523
1524
1525
1526//======================================================================
1527// Return an order in which to read the history such that _history[order[i]]
1528// will always correspond to the same set of consituent particles if
1529// two branching histories are equivalent in terms of the particles
1530// contained in any given pseudojet.
1531vector<int> ClusterSequence::unique_history_order() const {
1532
1533 // first construct an array that will tell us the lowest constituent
1534 // of a given jet -- this will always be one of the original
1535 // particles, whose order is well defined and so will help us to
1536 // follow the tree in a unique manner.
1537 valarray<int> lowest_constituent(_history.size());
1538 int hist_n = _history.size();
1539 lowest_constituent = hist_n; // give it a large number
1540 for (int i = 0; i < hist_n; i++) {
1541 // sets things up for the initial partons
1542 lowest_constituent[i] = min(lowest_constituent[i],i);
1543 // propagates them through to the children of this parton
1544 if (_history[i].child > 0) lowest_constituent[_history[i].child]
1545 = min(lowest_constituent[_history[i].child],lowest_constituent[i]);
1546 }
1547
1548 // establish an array for what we have and have not extracted so far
1549 valarray<bool> extracted(_history.size()); extracted = false;
1550 vector<int> unique_tree;
1551 unique_tree.reserve(_history.size());
1552
1553 // now work our way through the tree
1554 for (unsigned i = 0; i < n_particles(); i++) {
1555 if (!extracted[i]) {
1556 unique_tree.push_back(i);
1557 extracted[i] = true;
1558 _extract_tree_children(i, extracted, lowest_constituent, unique_tree);
1559 }
1560 }
1561
1562 return unique_tree;
1563}
1564
1565//======================================================================
1566// helper for unique_history_order
1567void ClusterSequence::_extract_tree_children(
1568 int position,
1569 valarray<bool> & extracted,
1570 const valarray<int> & lowest_constituent,
1571 vector<int> & unique_tree) const {
1572 if (!extracted[position]) {
1573 // that means we may have unidentified parents around, so go and
1574 // collect them (extracted[position]) will then be made true)
1575 _extract_tree_parents(position,extracted,lowest_constituent,unique_tree);
1576 }
1577
1578 // now look after the children...
1579 int child = _history[position].child;
1580 if (child >= 0) _extract_tree_children(child,extracted,lowest_constituent,unique_tree);
1581}
1582
1583
1584//======================================================================
1585// return the list of unclustered particles
1586vector<PseudoJet> ClusterSequence::unclustered_particles() const {
1587 vector<PseudoJet> unclustered;
1588 for (unsigned i = 0; i < n_particles() ; i++) {
1589 if (_history[i].child == Invalid)
1590 unclustered.push_back(_jets[_history[i].jetp_index]);
1591 }
1592 return unclustered;
1593}
1594
1595//======================================================================
1596/// Return the list of pseudojets in the ClusterSequence that do not
1597/// have children (and are not among the inclusive jets). They may
1598/// result from a clustering step or may be one of the pseudojets
1599/// returned by unclustered_particles().
1600vector<PseudoJet> ClusterSequence::childless_pseudojets() const {
1601 vector<PseudoJet> unclustered;
1602 for (unsigned i = 0; i < _history.size() ; i++) {
1603 if ((_history[i].child == Invalid) && (_history[i].parent2 != BeamJet))
1604 unclustered.push_back(_jets[_history[i].jetp_index]);
1605 }
1606 return unclustered;
1607}
1608
1609
1610
1611//----------------------------------------------------------------------
1612// returns true if the cluster sequence contains this jet (i.e. jet's
1613// structure is this cluster sequence's and the cluster history index
1614// is in a consistent range)
1615bool ClusterSequence::contains(const PseudoJet & jet) const {
1616 return jet.cluster_hist_index() >= 0
1617 && jet.cluster_hist_index() < int(_history.size())
1618 && jet.has_valid_cluster_sequence()
1619 && jet.associated_cluster_sequence() == this;
1620}
1621
1622
1623
1624//======================================================================
1625// helper for unique_history_order
1626void ClusterSequence::_extract_tree_parents(
1627 int position,
1628 valarray<bool> & extracted,
1629 const valarray<int> & lowest_constituent,
1630 vector<int> & unique_tree) const {
1631
1632 if (!extracted[position]) {
1633 int parent1 = _history[position].parent1;
1634 int parent2 = _history[position].parent2;
1635 // where relevant order parents so that we will first treat the
1636 // one containing the smaller "lowest_constituent"
1637 if (parent1 >= 0 && parent2 >= 0) {
1638 if (lowest_constituent[parent1] > lowest_constituent[parent2])
1639 std::swap(parent1, parent2);
1640 }
1641 // then actually run through the parents to extract the constituents...
1642 if (parent1 >= 0 && !extracted[parent1])
1643 _extract_tree_parents(parent1,extracted,lowest_constituent,unique_tree);
1644 if (parent2 >= 0 && !extracted[parent2])
1645 _extract_tree_parents(parent2,extracted,lowest_constituent,unique_tree);
1646 // finally declare this position to be accounted for and push it
1647 // onto our list.
1648 unique_tree.push_back(position);
1649 extracted[position] = true;
1650 }
1651}
1652
1653
1654//======================================================================
1655/// carries out the bookkeeping associated with the step of recombining
1656/// jet_i and jet_j (assuming a distance dij) and returns the index
1657/// of the recombined jet, newjet_k.
1658void ClusterSequence::_do_ij_recombination_step(
1659 const int jet_i, const int jet_j,
1660 const double dij,
1661 int & newjet_k) {
1662
1663 // Create the new jet by recombining the first two.
1664 //
1665 // For efficiency reasons, use a ctr that initialises only the
1666 // shared pointers, since the rest of the info will anyway be dealt
1667 // with by the recombiner.
1668 PseudoJet newjet(false);
1669 _jet_def.recombiner()->recombine(_jets[jet_i], _jets[jet_j], newjet);
1670 _jets.push_back(newjet);
1671 // original version...
1672 //_jets.push_back(_jets[jet_i] + _jets[jet_j]);
1673
1674 // get its index
1675 newjet_k = _jets.size()-1;
1676
1677 // get history index
1678 int newstep_k = _history.size();
1679 // and provide jet with the info
1680 _jets[newjet_k].set_cluster_hist_index(newstep_k);
1681
1682 // finally sort out the history
1683 int hist_i = _jets[jet_i].cluster_hist_index();
1684 int hist_j = _jets[jet_j].cluster_hist_index();
1685
1686 _add_step_to_history(min(hist_i, hist_j), max(hist_i,hist_j),
1687 newjet_k, dij);
1688
1689 // _add_step_to_history(newstep_k, min(hist_i, hist_j), max(hist_i,hist_j),
1690 // newjet_k, dij);
1691
1692
1693}
1694
1695
1696//======================================================================
1697/// carries out the bookkeeping associated with the step of recombining
1698/// jet_i with the beam
1699void ClusterSequence::_do_iB_recombination_step(
1700 const int jet_i, const double diB) {
1701 // recombine the jet with the beam
1702 _add_step_to_history(_jets[jet_i].cluster_hist_index(),BeamJet,
1703 Invalid, diB);
1704
1705 // // get history index
1706 // int newstep_k = _history.size();
1707 //
1708 // _add_step_to_history(newstep_k,_jets[jet_i].cluster_hist_index(),BeamJet,
1709 // Invalid, diB);
1710
1711}
1712
1713
1714// make sure the static member _changed_strategy_warning is defined.
1715LimitedWarning ClusterSequence::_changed_strategy_warning;
1716
1717
1718//----------------------------------------------------------------------
1719void ClusterSequence::_set_structure_shared_ptr(PseudoJet & j) {
1720 j.set_structure_shared_ptr(_structure_shared_ptr);
1721 // record the use count of the structure shared point to help
1722 // in case we want to ask the CS to handle its own memory
1723 _update_structure_use_count();
1724}
1725
1726
1727//----------------------------------------------------------------------
1728void ClusterSequence::_update_structure_use_count() {
1729 // record the use count of the structure shared point to help
1730 // in case we want to ask the CS to handle its own memory
1731 _structure_use_count_after_construction = _structure_shared_ptr.use_count();
1732}
1733
1734//----------------------------------------------------------------------
1735/// by calling this routine you tell the ClusterSequence to delete
1736/// itself when all the Pseudojets associated with it have gone out
1737/// of scope.
1738void ClusterSequence::delete_self_when_unused() {
1739 // the trick we use to handle this is to modify the use count;
1740 // that way the structure will be deleted when there are no external
1741 // objects left associated the CS and the structure's destructor will then
1742 // look after deleting the cluster sequence
1743
1744 // first make sure that there is at least one other object
1745 // associated with the CS
1746 int new_count = _structure_shared_ptr.use_count() - _structure_use_count_after_construction;
1747 if (new_count <= 0) {
1748 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");
1749 }
1750
1751 _structure_shared_ptr.set_count(new_count);
1752 _deletes_self_when_unused = true;
1753}
1754
1755
1756FASTJET_END_NAMESPACE
1757
Note: See TracBrowser for help on using the repository browser.