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