Fork me on GitHub

source: svn/trunk/Utilities/Fastjet/src/ClusterSequence.cc@ 217

Last change on this file since 217 was 213, checked in by severine ovyn, 16 years ago

remove banner

File size: 31.2 KB
Line 
1//STARTHEADER
2// $Id: ClusterSequence.cc,v 1.2 2009-01-29 23:28:55 ovyn Exp $
3//
4// Copyright (c) 2005-2006, Matteo Cacciari and Gavin Salam
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, write to the Free Software
26// Foundation, Inc.:
27// 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
28//----------------------------------------------------------------------
29//ENDHEADER
30
31#include "../include/fastjet/Error.hh"
32#include "../include/fastjet/PseudoJet.hh"
33#include "../include/fastjet/ClusterSequence.hh"
34#include "../include/fastjet/version.hh" // stores the current version number
35#include<iostream>
36#include<sstream>
37#include<cmath>
38#include<cstdlib>
39#include<cassert>
40#include<string>
41
42FASTJET_BEGIN_NAMESPACE // defined in fastjet/internal/base.hh
43
44using namespace std;
45
46//// initialised static member has to go in the .cc code
47JetAlgorithm ClusterSequence::_default_jet_algorithm = kt_algorithm;
48//
49
50
51//----------------------------------------------------------------------
52void ClusterSequence::_initialise_and_run (
53 const double & R,
54 const Strategy & strategy,
55 const bool & writeout_combinations) {
56
57 JetDefinition jet_def(_default_jet_algorithm, R, strategy);
58 _initialise_and_run(jet_def, writeout_combinations);
59}
60
61
62//----------------------------------------------------------------------
63void ClusterSequence::_initialise_and_run (
64 const JetDefinition & jet_def,
65 const bool & writeout_combinations) {
66
67 // transfer all relevant info into internal variables
68 _decant_options(jet_def, writeout_combinations);
69
70 // set up the history entries for the initial particles (those
71 // currently in _jets)
72 _fill_initial_history();
73
74 // run the plugin if that's what's decreed
75 if (_jet_algorithm == plugin_algorithm) {
76 _plugin_activated = true;
77 _jet_def.plugin()->run_clustering( (*this) );
78 _plugin_activated = false;
79 return;
80 }
81
82
83 // automatically redefine the strategy according to N if that is
84 // what the user requested -- transition points (and especially
85 // their R-dependence) are based on empirical observations for a
86 // R=0.4, 0.7 and 1.0, running on toth (3.4GHz, Pentium IV D [dual
87 // core] with 2MB of cache).
88 if (_strategy == Best) {
89 int N = _jets.size();
90 if (N > 6200/pow(_Rparam,2.0)
91 && jet_def.jet_algorithm() == cambridge_algorithm) {
92 _strategy = NlnNCam;}
93 else
94#ifndef DROP_CGAL
95 if (N > 16000/pow(_Rparam,1.15)) {
96 _strategy = NlnN; }
97 else
98#endif // DROP_CGAL
99 if (N > 450) {
100 _strategy = N2MinHeapTiled;
101 }
102 else if (N > 55*max(0.5,min(1.0,_Rparam))) {// empirical scaling with R
103 _strategy = N2Tiled;
104 } else {
105 _strategy = N2Plain;
106 }
107 }
108
109
110 // run the code containing the selected strategy
111 if (_strategy == NlnN || _strategy == NlnN3pi
112 || _strategy == NlnN4pi ) {
113 this->_delaunay_cluster();
114 } else if (_strategy == N3Dumb ) {
115 this->_really_dumb_cluster();
116 } else if (_strategy == N2Tiled) {
117 this->_faster_tiled_N2_cluster();
118 } else if (_strategy == N2PoorTiled) {
119 this->_tiled_N2_cluster();
120 } else if (_strategy == N2Plain) {
121 this->_simple_N2_cluster();
122 } else if (_strategy == N2MinHeapTiled) {
123 this->_minheap_faster_tiled_N2_cluster();
124 } else if (_strategy == NlnNCam4pi) {
125 this->_CP2DChan_cluster();
126 } else if (_strategy == NlnNCam2pi2R) {
127 this->_CP2DChan_cluster_2pi2R();
128 } else if (_strategy == NlnNCam) {
129 this->_CP2DChan_cluster_2piMultD();
130 } else {
131 ostringstream err;
132 err << "Unrecognised value for strategy: "<<_strategy;
133 throw Error(err.str());
134 //assert(false);
135 }
136}
137
138
139// these needs to be defined outside the class definition.
140bool ClusterSequence::_first_time = true;
141int ClusterSequence::_n_exclusive_warnings = 0;
142
143
144//----------------------------------------------------------------------
145// the version string
146string fastjet_version_string() {
147 return "FastJet version "+string(fastjet_version);
148}
149
150
151//----------------------------------------------------------------------
152// prints a banner on the first call
153void ClusterSequence::_print_banner() {
154
155 if (!_first_time) {return;}
156 _first_time = false;
157
158
159 //Symp. Discr. Alg, p.472 (2002) and CGAL (http://www.cgal.org);
160
161 cout << "#--------------------------------------------------------------------------\n";
162 cout << "# FastJet release " << fastjet_version << endl;
163 cout << "# Written by M. Cacciari, G.P. Salam and G. Soyez \n";
164 cout << "# http://www.lpthe.jussieu.fr/~salam/fastjet \n";
165 cout << "# \n";
166 cout << "# Longitudinally invariant Kt, anti-Kt, and inclusive Cambridge/Aachen \n";
167 cout << "# clustering using fast geometric algorithms, with area measures and optional\n";
168 cout << "# external jet-finder plugins. \n";
169 cout << "# Please cite Phys. Lett. B641 (2006) [hep-ph/0512210] if you use this code.\n";
170 cout << "# \n";
171 cout << "# This package uses T.Chan's closest pair algorithm, Proc.13th ACM-SIAM \n";
172 cout << "# Symp. Discr. Alg, p.472 (2002), S.Fortune's Voronoi algorithm and code " ;
173#ifndef DROP_CGAL
174 cout << endl << "# and CGAL: http://www.cgal.org/";
175#endif // DROP_CGAL
176 cout << ".\n";
177 cout << "#-------------------------------------------------------------------------\n";
178}
179
180//----------------------------------------------------------------------
181// transfer all relevant info into internal variables
182void ClusterSequence::_decant_options(const JetDefinition & jet_def,
183 const bool & writeout_combinations) {
184
185 // let the user know what's going on
186// _print_banner();
187
188 // make a local copy of the jet definition (for future use?)
189 _jet_def = jet_def;
190
191 _writeout_combinations = writeout_combinations;
192 _jet_algorithm = jet_def.jet_algorithm();
193 _Rparam = jet_def.R(); _R2 = _Rparam*_Rparam; _invR2 = 1.0/_R2;
194 _strategy = jet_def.strategy();
195
196 // disallow interference from the plugin
197 _plugin_activated = false;
198
199}
200
201
202//----------------------------------------------------------------------
203// initialise the history in a standard way
204void ClusterSequence::_fill_initial_history () {
205
206 if (_jets.size() == 0) {throw Error("Cannot run jet-finder on empty event");}
207
208 // reserve sufficient space for everything
209 _jets.reserve(_jets.size()*2);
210 _history.reserve(_jets.size()*2);
211
212 for (int i = 0; i < static_cast<int>(_jets.size()) ; i++) {
213 history_element element;
214 element.parent1 = InexistentParent;
215 element.parent2 = InexistentParent;
216 element.child = Invalid;
217 element.jetp_index = i;
218 element.dij = 0.0;
219 element.max_dij_so_far = 0.0;
220
221 _history.push_back(element);
222
223 // do any momentum preprocessing needed by the recombination scheme
224 _jet_def.recombiner()->preprocess(_jets[i]);
225
226 // get cross-referencing right from PseudoJets
227 _jets[i].set_cluster_hist_index(i);
228 }
229 _initial_n = _jets.size();
230}
231
232
233//----------------------------------------------------------------------
234// Return the component corresponding to the specified index.
235// taken from CLHEP
236string ClusterSequence::strategy_string () const {
237 string strategy;
238 switch(_strategy) {
239 case NlnN:
240 strategy = "NlnN"; break;
241 case NlnN3pi:
242 strategy = "NlnN3pi"; break;
243 case NlnN4pi:
244 strategy = "NlnN4pi"; break;
245 case N2Plain:
246 strategy = "N2Plain"; break;
247 case N2Tiled:
248 strategy = "N2Tiled"; break;
249 case N2MinHeapTiled:
250 strategy = "N2MinHeapTiled"; break;
251 case N2PoorTiled:
252 strategy = "N2PoorTiled"; break;
253 case N3Dumb:
254 strategy = "N3Dumb"; break;
255 case NlnNCam4pi:
256 strategy = "NlnNCam4pi"; break;
257 case NlnNCam2pi2R:
258 strategy = "NlnNCam2pi2R"; break;
259 case NlnNCam:
260 strategy = "NlnNCam"; break; // 2piMultD
261 case plugin_strategy:
262 strategy = "plugin strategy"; break;
263 default:
264 strategy = "Unrecognized";
265 }
266 return strategy;
267}
268
269
270double ClusterSequence::jet_scale_for_algorithm(
271 const PseudoJet & jet) const {
272 if (_jet_algorithm == kt_algorithm) {return jet.kt2();}
273 else if (_jet_algorithm == cambridge_algorithm) {return 1.0;}
274 else if (_jet_algorithm == antikt_algorithm) {
275 double kt2=jet.kt2();
276 return kt2 > 1e-300 ? 1.0/kt2 : 1e300;
277 } else if (_jet_algorithm == cambridge_for_passive_algorithm) {
278 double kt2 = jet.kt2();
279 double lim = _jet_def.extra_param();
280 if (kt2 < lim*lim && kt2 != 0.0) {
281 return 1.0/kt2;
282 } else {return 1.0;}
283 } else {throw Error("Unrecognised jet finder");}
284}
285
286
287//----------------------------------------------------------------------
288/// transfer the sequence contained in other_seq into our own;
289/// any plugin "extras" contained in the from_seq will be lost
290/// from there.
291void ClusterSequence::transfer_from_sequence(ClusterSequence & from_seq) {
292
293 // the metadata
294 _jet_def = from_seq._jet_def ;
295 _writeout_combinations = from_seq._writeout_combinations ;
296 _initial_n = from_seq._initial_n ;
297 _Rparam = from_seq._Rparam ;
298 _R2 = from_seq._R2 ;
299 _invR2 = from_seq._invR2 ;
300 _strategy = from_seq._strategy ;
301 _jet_algorithm = from_seq._jet_algorithm ;
302 _plugin_activated = from_seq._plugin_activated ;
303
304 // the data
305 _jets = from_seq._jets;
306 _history = from_seq._history;
307 // the following transferse ownership of the extras from the from_seq
308 _extras = from_seq._extras;
309
310}
311
312//----------------------------------------------------------------------
313// record an ij recombination and reset the _jets[newjet_k] momentum and
314// user index to be those of newjet
315void ClusterSequence::plugin_record_ij_recombination(
316 int jet_i, int jet_j, double dij,
317 const PseudoJet & newjet, int & newjet_k) {
318
319 plugin_record_ij_recombination(jet_i, jet_j, dij, newjet_k);
320
321 // now transfer newjet into place
322 int tmp_index = _jets[newjet_k].cluster_hist_index();
323 _jets[newjet_k] = newjet;
324 _jets[newjet_k].set_cluster_hist_index(tmp_index);
325}
326
327
328//----------------------------------------------------------------------
329// return all inclusive jets with pt > ptmin
330vector<PseudoJet> ClusterSequence::inclusive_jets (const double & ptmin) const{
331 double dcut = ptmin*ptmin;
332 int i = _history.size() - 1; // last jet
333 vector<PseudoJet> jets;
334 if (_jet_algorithm == kt_algorithm) {
335 while (i >= 0) {
336 // with our specific definition of dij and diB (i.e. R appears only in
337 // dij), then dij==diB is the same as the jet.perp2() and we can exploit
338 // this in selecting the jets...
339 if (_history[i].max_dij_so_far < dcut) {break;}
340 if (_history[i].parent2 == BeamJet && _history[i].dij >= dcut) {
341 // for beam jets
342 int parent1 = _history[i].parent1;
343 jets.push_back(_jets[_history[parent1].jetp_index]);}
344 i--;
345 }
346 } else if (_jet_algorithm == cambridge_algorithm) {
347 while (i >= 0) {
348 // inclusive jets are all at end of clustering sequence in the
349 // Cambridge algorithm -- so if we find a non-exclusive jet, then
350 // we can exit
351 if (_history[i].parent2 != BeamJet) {break;}
352 int parent1 = _history[i].parent1;
353 const PseudoJet & jet = _jets[_history[parent1].jetp_index];
354 if (jet.perp2() >= dcut) {jets.push_back(jet);}
355 i--;
356 }
357 } else if (_jet_algorithm == plugin_algorithm
358 || _jet_algorithm == antikt_algorithm
359 || _jet_algorithm == cambridge_for_passive_algorithm) {
360 // for inclusive jets with a plugin algorithm, we make no
361 // assumptions about anything (relation of dij to momenta,
362 // ordering of the dij, etc.)
363 while (i >= 0) {
364 if (_history[i].parent2 == BeamJet) {
365 int parent1 = _history[i].parent1;
366 const PseudoJet & jet = _jets[_history[parent1].jetp_index];
367 if (jet.perp2() >= dcut) {jets.push_back(jet);}
368 }
369 i--;
370 }
371 } else {throw Error("Unrecognized jet algorithm");}
372 return jets;
373}
374
375
376//----------------------------------------------------------------------
377// return the number of exclusive jets that would have been obtained
378// running the algorithm in exclusive mode with the given dcut
379int ClusterSequence::n_exclusive_jets (const double & dcut) const {
380
381 // first locate the point where clustering would have stopped (i.e. the
382 // first time max_dij_so_far > dcut)
383 int i = _history.size() - 1; // last jet
384 while (i >= 0) {
385 if (_history[i].max_dij_so_far <= dcut) {break;}
386 i--;
387 }
388 int stop_point = i + 1;
389 // relation between stop_point, njets assumes one extra jet disappears
390 // at each clustering.
391 int njets = 2*_initial_n - stop_point;
392 return njets;
393}
394
395//----------------------------------------------------------------------
396// return all exclusive jets that would have been obtained running
397// the algorithm in exclusive mode with the given dcut
398vector<PseudoJet> ClusterSequence::exclusive_jets (const double & dcut) const {
399 int njets = n_exclusive_jets(dcut);
400 return exclusive_jets(njets);
401}
402
403
404//----------------------------------------------------------------------
405// return the jets obtained by clustering the event to n jets.
406vector<PseudoJet> ClusterSequence::exclusive_jets (const int & njets) const {
407
408 // make sure the user does not ask for more than jets than there
409 // were particles in the first place.
410 assert (njets <= _initial_n);
411
412 // provide a warning when extracting exclusive jets
413 if (_jet_def.jet_algorithm() != kt_algorithm && _n_exclusive_warnings < 5) {
414 _n_exclusive_warnings++;
415 cerr << "FastJet WARNING: dcut and exclusive jets for jet-finders other than kt should be interpreted with care." << endl;
416 }
417
418
419 // calculate the point where we have to stop the clustering.
420 // relation between stop_point, njets assumes one extra jet disappears
421 // at each clustering.
422 int stop_point = 2*_initial_n - njets;
423
424 // some sanity checking to make sure that e+e- does not give us
425 // surprises (should we ever implement e+e-)...
426 if (2*_initial_n != static_cast<int>(_history.size())) {
427 ostringstream err;
428 err << "2*_initial_n != _history.size() -- this endangers internal assumptions!\n";
429 throw Error(err.str());
430 //assert(false);
431 }
432
433 // now go forwards and reconstitute the jets that we have --
434 // basically for any history element, see if the parent jets to
435 // which it refers were created before the stopping point -- if they
436 // were then add them to the list, otherwise they are subsequent
437 // recombinations of the jets that we are looking for.
438 vector<PseudoJet> jets;
439 for (unsigned int i = stop_point; i < _history.size(); i++) {
440 int parent1 = _history[i].parent1;
441 if (parent1 < stop_point) {
442 jets.push_back(_jets[_history[parent1].jetp_index]);
443 }
444 int parent2 = _history[i].parent2;
445 if (parent2 < stop_point && parent2 > 0) {
446 jets.push_back(_jets[_history[parent2].jetp_index]);
447 }
448
449 }
450
451 // sanity check...
452 if (static_cast<int>(jets.size()) != njets) {
453 ostringstream err;
454 err << "ClusterSequence::exclusive_jets: size of returned vector ("
455 <<jets.size()<<") does not coincide with requested number of jets ("
456 <<njets<<")";
457 throw Error(err.str());
458 }
459
460 return jets;
461}
462
463//----------------------------------------------------------------------
464/// return the dmin corresponding to the recombination that went from
465/// n+1 to n jets
466double ClusterSequence::exclusive_dmerge (const int & njets) const {
467 assert(njets >= 0);
468 if (njets >= _initial_n) {return 0.0;}
469 return _history[2*_initial_n-njets-1].dij;
470}
471
472
473//----------------------------------------------------------------------
474/// return the maximum of the dmin encountered during all recombinations
475/// up to the one that led to an n-jet final state; identical to
476/// exclusive_dmerge, except in cases where the dmin do not increase
477/// monotonically.
478double ClusterSequence::exclusive_dmerge_max (const int & njets) const {
479 assert(njets >= 0);
480 if (njets >= _initial_n) {return 0.0;}
481 return _history[2*_initial_n-njets-1].max_dij_so_far;
482}
483
484
485//----------------------------------------------------------------------
486// work through the object's history until
487bool ClusterSequence::object_in_jet(const PseudoJet & object,
488 const PseudoJet & jet) const {
489
490 // make sure the object conceivably belongs to this clustering
491 // sequence
492 assert(_potentially_valid(object) && _potentially_valid(jet));
493
494 const PseudoJet * this_object = &object;
495 const PseudoJet * childp;
496 while(true) {
497 if (this_object->cluster_hist_index() == jet.cluster_hist_index()) {
498 return true;
499 } else if (has_child(*this_object, childp)) {this_object = childp;}
500 else {return false;}
501 }
502}
503
504//----------------------------------------------------------------------
505/// if the jet has parents in the clustering, it returns true
506/// and sets parent1 and parent2 equal to them.
507///
508/// if it has no parents it returns false and sets parent1 and
509/// parent2 to zero
510bool ClusterSequence::has_parents(const PseudoJet & jet, PseudoJet & parent1,
511 PseudoJet & parent2) const {
512
513 const history_element & hist = _history[jet.cluster_hist_index()];
514
515 // make sure we do not run into any unexpected situations --
516 // i.e. both parents valid, or neither
517 assert ((hist.parent1 >= 0 && hist.parent2 >= 0) ||
518 (hist.parent1 < 0 && hist.parent2 < 0));
519
520 if (hist.parent1 < 0) {
521 parent1 = PseudoJet(0.0,0.0,0.0,0.0);
522 parent2 = parent1;
523 return false;
524 } else {
525 parent1 = _jets[_history[hist.parent1].jetp_index];
526 parent2 = _jets[_history[hist.parent2].jetp_index];
527 // order the parents in decreasing pt
528 if (parent1.perp2() < parent2.perp2()) swap(parent1,parent2);
529 return true;
530 }
531}
532
533//----------------------------------------------------------------------
534/// if the jet has a child then return true and give the child jet
535/// otherwise return false and set the child to zero
536bool ClusterSequence::has_child(const PseudoJet & jet, PseudoJet & child) const {
537
538 //const history_element & hist = _history[jet.cluster_hist_index()];
539 //
540 //if (hist.child >= 0) {
541 // child = _jets[_history[hist.child].jetp_index];
542 // return true;
543 //} else {
544 // child = PseudoJet(0.0,0.0,0.0,0.0);
545 // return false;
546 //}
547 const PseudoJet * childp;
548 bool res = has_child(jet, childp);
549 if (res) {
550 child = *childp;
551 return true;
552 } else {
553 child = PseudoJet(0.0,0.0,0.0,0.0);
554 return false;
555 }
556}
557
558bool ClusterSequence::has_child(const PseudoJet & jet, const PseudoJet * & childp) const {
559
560 const history_element & hist = _history[jet.cluster_hist_index()];
561
562 // check that this jet has a child and that the child corresponds to
563 // a true jet [RETHINK-IF-CHANGE-NUMBERING: what is the right
564 // behaviour if the child is the same jet but made inclusive...?]
565 if (hist.child >= 0 && _history[hist.child].jetp_index >= 0) {
566 childp = &(_jets[_history[hist.child].jetp_index]);
567 return true;
568 } else {
569 childp = NULL;
570 return false;
571 }
572}
573
574
575//----------------------------------------------------------------------
576/// if this jet has a child (and so a partner) return true
577/// and give the partner, otherwise return false and set the
578/// partner to zero
579bool ClusterSequence::has_partner(const PseudoJet & jet,
580 PseudoJet & partner) const {
581
582 const history_element & hist = _history[jet.cluster_hist_index()];
583
584 // make sure we have a child and that the child does not correspond
585 // to a clustering with the beam (or some other invalid quantity)
586 if (hist.child >= 0 && _history[hist.child].parent2 >= 0) {
587 const history_element & child_hist = _history[hist.child];
588 if (child_hist.parent1 == jet.cluster_hist_index()) {
589 // partner will be child's parent2 -- for iB clustering
590 // parent2 will not be valid
591 partner = _jets[_history[child_hist.parent2].jetp_index];
592 } else {
593 // partner will be child's parent1
594 partner = _jets[_history[child_hist.parent1].jetp_index];
595 }
596 return true;
597 } else {
598 partner = PseudoJet(0.0,0.0,0.0,0.0);
599 return false;
600 }
601}
602
603
604//----------------------------------------------------------------------
605// return a vector of the particles that make up a jet
606vector<PseudoJet> ClusterSequence::constituents (const PseudoJet & jet) const {
607 vector<PseudoJet> subjets;
608 add_constituents(jet, subjets);
609 return subjets;
610}
611
612//----------------------------------------------------------------------
613/// output the supplied vector of jets in a format that can be read
614/// by an appropriate root script; the format is:
615/// jet-n jet-px jet-py jet-pz jet-E
616/// particle-n particle-rap particle-phi particle-pt
617/// particle-n particle-rap particle-phi particle-pt
618/// ...
619/// #END
620/// ... [i.e. above repeated]
621void ClusterSequence::print_jets_for_root(const std::vector<PseudoJet> & jets,
622 ostream & ostr) const {
623 for (unsigned i = 0; i < jets.size(); i++) {
624 ostr << i << " "
625 << jets[i].px() << " "
626 << jets[i].py() << " "
627 << jets[i].pz() << " "
628 << jets[i].E() << endl;
629 vector<PseudoJet> cst = constituents(jets[i]);
630 for (unsigned j = 0; j < cst.size() ; j++) {
631 ostr << " " << j << " "
632 << cst[j].rap() << " "
633 << cst[j].phi() << " "
634 << cst[j].perp() << endl;
635 }
636 ostr << "#END" << endl;
637 }
638}
639
640// Not yet. Perhaps in a future release
641// //----------------------------------------------------------------------
642// // print out all inclusive jets with pt > ptmin
643// void ClusterSequence::print_jets (const double & ptmin) const{
644// vector<PseudoJet> jets = sorted_by_pt(inclusive_jets(ptmin));
645//
646// for (size_t j = 0; j < jets.size(); j++) {
647// printf("%5u %7.3f %7.3f %9.3f\n",
648// j,jets[j].rap(),jets[j].phi(),jets[j].perp());
649// }
650// }
651
652//----------------------------------------------------------------------
653/// returns a vector of size n_particles() which indicates, for
654/// each of the initial particles (in the order in which they were
655/// supplied), which of the supplied jets it belongs to; if it does
656/// not belong to any of the supplied jets, the index is set to -1;
657vector<int> ClusterSequence::particle_jet_indices(
658 const vector<PseudoJet> & jets) const {
659
660 vector<int> indices(n_particles());
661
662 // first label all particles as not belonging to any jets
663 for (unsigned ipart = 0; ipart < n_particles(); ipart++)
664 indices[ipart] = -1;
665
666 // then for each of the jets relabel its consituents as belonging to
667 // that jet
668 for (unsigned ijet = 0; ijet < jets.size(); ijet++) {
669
670 vector<PseudoJet> jet_constituents(constituents(jets[ijet]));
671
672 for (unsigned ip = 0; ip < jet_constituents.size(); ip++) {
673 // a safe (if slightly redundant) way of getting the particle
674 // index (for initial particles it is actually safe to assume
675 // ipart=iclust).
676 unsigned iclust = jet_constituents[ip].cluster_hist_index();
677 unsigned ipart = history()[iclust].jetp_index;
678 indices[ipart] = ijet;
679 }
680 }
681
682 return indices;
683}
684
685
686//----------------------------------------------------------------------
687// recursive routine that adds on constituents of jet to the subjet_vector
688void ClusterSequence::add_constituents (
689 const PseudoJet & jet, vector<PseudoJet> & subjet_vector) const {
690 // find out position in cluster history
691 int i = jet.cluster_hist_index();
692 int parent1 = _history[i].parent1;
693 int parent2 = _history[i].parent2;
694
695 if (parent1 == InexistentParent) {
696 // It is an original particle (labelled by its parent having value
697 // InexistentParent), therefore add it on to the subjet vector
698 subjet_vector.push_back(jet);
699 return;
700 }
701
702 // add parent 1
703 add_constituents(_jets[_history[parent1].jetp_index], subjet_vector);
704
705 // see if parent2 is a real jet; if it is then add its constituents
706 if (parent2 != BeamJet) {
707 add_constituents(_jets[_history[parent2].jetp_index], subjet_vector);
708 }
709}
710
711
712
713//----------------------------------------------------------------------
714// initialise the history in a standard way
715void ClusterSequence::_add_step_to_history (
716 const int & step_number, const int & parent1,
717 const int & parent2, const int & jetp_index,
718 const double & dij) {
719
720 history_element element;
721 element.parent1 = parent1;
722 element.parent2 = parent2;
723 element.jetp_index = jetp_index;
724 element.child = Invalid;
725 element.dij = dij;
726 element.max_dij_so_far = max(dij,_history[_history.size()-1].max_dij_so_far);
727 _history.push_back(element);
728
729 int local_step = _history.size()-1;
730 assert(local_step == step_number);
731
732 assert(parent1 >= 0);
733 _history[parent1].child = local_step;
734 if (parent2 >= 0) {_history[parent2].child = local_step;}
735
736 // get cross-referencing right from PseudoJets
737 if (jetp_index != Invalid) {
738 assert(jetp_index >= 0);
739 //cout << _jets.size() <<" "<<jetp_index<<"\n";
740 _jets[jetp_index].set_cluster_hist_index(local_step);
741 }
742
743 if (_writeout_combinations) {
744 cout << local_step << ": "
745 << parent1 << " with " << parent2
746 << "; y = "<< dij<<endl;
747 }
748
749}
750
751
752
753
754//======================================================================
755// Return an order in which to read the history such that _history[order[i]]
756// will always correspond to the same set of consituent particles if
757// two branching histories are equivalent in terms of the particles
758// contained in any given pseudojet.
759vector<int> ClusterSequence::unique_history_order() const {
760
761 // first construct an array that will tell us the lowest constituent
762 // of a given jet -- this will always be one of the original
763 // particles, whose order is well defined and so will help us to
764 // follow the tree in a unique manner.
765 valarray<int> lowest_constituent(_history.size());
766 int hist_n = _history.size();
767 lowest_constituent = hist_n; // give it a large number
768 for (int i = 0; i < hist_n; i++) {
769 // sets things up for the initial partons
770 lowest_constituent[i] = min(lowest_constituent[i],i);
771 // propagates them through to the children of this parton
772 if (_history[i].child > 0) lowest_constituent[_history[i].child]
773 = min(lowest_constituent[_history[i].child],lowest_constituent[i]);
774 }
775
776 // establish an array for what we have and have not extracted so far
777 valarray<bool> extracted(_history.size()); extracted = false;
778 vector<int> unique_tree;
779 unique_tree.reserve(_history.size());
780
781 // now work our way through the tree
782 for (unsigned i = 0; i < n_particles(); i++) {
783 if (!extracted[i]) {
784 unique_tree.push_back(i);
785 extracted[i] = true;
786 _extract_tree_children(i, extracted, lowest_constituent, unique_tree);
787 }
788 }
789
790 return unique_tree;
791}
792
793//======================================================================
794// helper for unique_history_order
795void ClusterSequence::_extract_tree_children(
796 int position,
797 valarray<bool> & extracted,
798 const valarray<int> & lowest_constituent,
799 vector<int> & unique_tree) const {
800 if (!extracted[position]) {
801 // that means we may have unidentified parents around, so go and
802 // collect them (extracted[position]) will then be made true)
803 _extract_tree_parents(position,extracted,lowest_constituent,unique_tree);
804 }
805
806 // now look after the children...
807 int child = _history[position].child;
808 if (child >= 0) _extract_tree_children(child,extracted,lowest_constituent,unique_tree);
809}
810
811
812//======================================================================
813// return the list of unclustered particles
814vector<PseudoJet> ClusterSequence::unclustered_particles() const {
815 vector<PseudoJet> unclustered;
816 for (unsigned i = 0; i < n_particles() ; i++) {
817 if (_history[i].child == Invalid)
818 unclustered.push_back(_jets[_history[i].jetp_index]);
819 }
820 return unclustered;
821}
822
823
824
825//======================================================================
826// helper for unique_history_order
827void ClusterSequence::_extract_tree_parents(
828 int position,
829 valarray<bool> & extracted,
830 const valarray<int> & lowest_constituent,
831 vector<int> & unique_tree) const {
832
833 if (!extracted[position]) {
834 int parent1 = _history[position].parent1;
835 int parent2 = _history[position].parent2;
836 // where relevant order parents so that we will first treat the
837 // one containing the smaller "lowest_constituent"
838 if (parent1 >= 0 && parent2 >= 0) {
839 if (lowest_constituent[parent1] > lowest_constituent[parent2])
840 swap(parent1, parent2);
841 }
842 // then actually run through the parents to extract the constituents...
843 if (parent1 >= 0 && !extracted[parent1])
844 _extract_tree_parents(parent1,extracted,lowest_constituent,unique_tree);
845 if (parent2 >= 0 && !extracted[parent2])
846 _extract_tree_parents(parent2,extracted,lowest_constituent,unique_tree);
847 // finally declare this position to be accounted for and push it
848 // onto our list.
849 unique_tree.push_back(position);
850 extracted[position] = true;
851 }
852}
853
854
855//======================================================================
856/// carries out the bookkeeping associated with the step of recombining
857/// jet_i and jet_j (assuming a distance dij) and returns the index
858/// of the recombined jet, newjet_k.
859void ClusterSequence::_do_ij_recombination_step(
860 const int & jet_i, const int & jet_j,
861 const double & dij,
862 int & newjet_k) {
863
864 // create the new jet by recombining the first two
865 PseudoJet newjet;
866 _jet_def.recombiner()->recombine(_jets[jet_i], _jets[jet_j], newjet);
867 _jets.push_back(newjet);
868 // original version...
869 //_jets.push_back(_jets[jet_i] + _jets[jet_j]);
870
871 // get its index
872 newjet_k = _jets.size()-1;
873
874 // get history index
875 int newstep_k = _history.size();
876 // and provide jet with the info
877 _jets[newjet_k].set_cluster_hist_index(newstep_k);
878
879 // finally sort out the history
880 int hist_i = _jets[jet_i].cluster_hist_index();
881 int hist_j = _jets[jet_j].cluster_hist_index();
882
883 _add_step_to_history(newstep_k, min(hist_i, hist_j), max(hist_i,hist_j),
884 newjet_k, dij);
885
886}
887
888
889//======================================================================
890/// carries out the bookkeeping associated with the step of recombining
891/// jet_i with the beam
892void ClusterSequence::_do_iB_recombination_step(
893 const int & jet_i, const double & diB) {
894 // get history index
895 int newstep_k = _history.size();
896
897 // recombine the jet with the beam
898 _add_step_to_history(newstep_k,_jets[jet_i].cluster_hist_index(),BeamJet,
899 Invalid, diB);
900
901}
902
903FASTJET_END_NAMESPACE
904
Note: See TracBrowser for help on using the repository browser.