Fork me on GitHub

source: git/external/fastjet/ClusterSequence.cc@ 4b9a2dc

ImprovedOutputFile Timing dual_readout llp
Last change on this file since 4b9a2dc was 35cdc46, checked in by Pavel Demin <demin@…>, 10 years ago

upgrade FastJet to version 3.1.0-beta.1, upgrade Nsubjettiness to version 2.1.0, add SoftKiller version 1.0.0

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