00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013
00014
00015
00016
00017
00018
00019
00020
00021
00022
00023
00024
00025
00026
00027
00028
00029
00030
00031 #include "fastjet/Error.hh"
00032 #include "fastjet/PseudoJet.hh"
00033 #include "fastjet/ClusterSequence.hh"
00034 #include "fastjet/version.hh"
00035 #include<iostream>
00036 #include<sstream>
00037 #include<fstream>
00038 #include<cmath>
00039 #include<cstdlib>
00040 #include<cassert>
00041 #include<string>
00042 #include<set>
00043
00044 FASTJET_BEGIN_NAMESPACE
00045
00046 using namespace std;
00047
00049 JetAlgorithm ClusterSequence::_default_jet_algorithm = kt_algorithm;
00050
00051
00052
00053
00054 ClusterSequence::~ClusterSequence () {}
00055
00056
00057 void ClusterSequence::_initialise_and_run (
00058 const double & R,
00059 const Strategy & strategy,
00060 const bool & writeout_combinations) {
00061
00062 JetDefinition jet_def(_default_jet_algorithm, R, strategy);
00063 _initialise_and_run(jet_def, writeout_combinations);
00064 }
00065
00066
00067
00068 void ClusterSequence::_initialise_and_run (
00069 const JetDefinition & jet_def,
00070 const bool & writeout_combinations) {
00071
00072
00073 _decant_options(jet_def, writeout_combinations);
00074
00075
00076
00077 _fill_initial_history();
00078
00079
00080 if (n_particles() == 0) return;
00081
00082
00083 if (_jet_algorithm == plugin_algorithm) {
00084
00085 _plugin_activated = true;
00086
00087 _jet_def.plugin()->run_clustering( (*this) );
00088 _plugin_activated = false;
00089 return;
00090 } else if (_jet_algorithm == ee_kt_algorithm ||
00091 _jet_algorithm == ee_genkt_algorithm) {
00092
00093 _strategy = N2Plain;
00094 if (_jet_algorithm == ee_kt_algorithm) {
00095
00096
00097
00098 assert(_Rparam > 2.0);
00099
00100
00101
00102 _invR2 = 1.0;
00103 } else {
00104
00105
00106
00107
00108 if (_Rparam > pi) {
00109
00110
00111
00112 _R2 = 2 * ( 3.0 + cos(_Rparam) );
00113 } else {
00114 _R2 = 2 * ( 1.0 - cos(_Rparam) );
00115 }
00116 _invR2 = 1.0/_R2;
00117 }
00118
00119 _simple_N2_cluster_EEBriefJet();
00120 return;
00121 }
00122
00123
00124
00125
00126
00127
00128
00129 if (_strategy == Best) {
00130 int N = _jets.size();
00131 if (N > 6200/pow(_Rparam,2.0)
00132 && jet_def.jet_algorithm() == cambridge_algorithm) {
00133 _strategy = NlnNCam;}
00134 else
00135 #ifndef DROP_CGAL
00136 if ((N > 16000/pow(_Rparam,1.15) && jet_def.jet_algorithm() != antikt_algorithm)
00137 || N > 35000/pow(_Rparam,1.15)) {
00138 _strategy = NlnN; }
00139 else
00140 #endif // DROP_CGAL
00141 if (N > 450) {
00142 _strategy = N2MinHeapTiled;
00143 }
00144 else if (N > 55*max(0.5,min(1.0,_Rparam))) {
00145 _strategy = N2Tiled;
00146 } else {
00147 _strategy = N2Plain;
00148 }
00149 }
00150
00151
00152
00153 if (_strategy == NlnN || _strategy == NlnN3pi
00154 || _strategy == NlnN4pi ) {
00155 this->_delaunay_cluster();
00156 } else if (_strategy == N3Dumb ) {
00157 this->_really_dumb_cluster();
00158 } else if (_strategy == N2Tiled) {
00159 this->_faster_tiled_N2_cluster();
00160 } else if (_strategy == N2PoorTiled) {
00161 this->_tiled_N2_cluster();
00162 } else if (_strategy == N2Plain) {
00163
00164
00165 this->_simple_N2_cluster_BriefJet();
00166 } else if (_strategy == N2MinHeapTiled) {
00167 this->_minheap_faster_tiled_N2_cluster();
00168 } else if (_strategy == NlnNCam4pi) {
00169 this->_CP2DChan_cluster();
00170 } else if (_strategy == NlnNCam2pi2R) {
00171 this->_CP2DChan_cluster_2pi2R();
00172 } else if (_strategy == NlnNCam) {
00173 this->_CP2DChan_cluster_2piMultD();
00174 } else {
00175 ostringstream err;
00176 err << "Unrecognised value for strategy: "<<_strategy;
00177 throw Error(err.str());
00178
00179 }
00180 }
00181
00182
00183
00184 bool ClusterSequence::_first_time = true;
00185 int ClusterSequence::_n_exclusive_warnings = 0;
00186
00187
00188
00189
00190 string fastjet_version_string() {
00191 return "FastJet version "+string(fastjet_version);
00192 }
00193
00194
00195
00196
00197 void ClusterSequence::_print_banner() {
00198
00199 if (!_first_time) {return;}
00200 _first_time = false;
00201
00202
00203
00204
00205 cout << "#--------------------------------------------------------------------------\n";
00206 cout << "# FastJet release " << fastjet_version << endl;
00207 cout << "# Written by M. Cacciari, G.P. Salam and G. Soyez \n";
00208 cout << "# http://www.fastjet.fr \n";
00209 cout << "# \n";
00210 cout << "# Longitudinally invariant Kt, anti-Kt, and inclusive Cambridge/Aachen \n";
00211 cout << "# clustering using fast geometric algorithms, with area measures and optional\n";
00212 cout << "# external jet-finder plugins. \n";
00213 cout << "# Please cite Phys. Lett. B641 (2006) [hep-ph/0512210] if you use this code.\n";
00214 cout << "# \n";
00215 cout << "# This package uses T.Chan's closest pair algorithm, Proc.13th ACM-SIAM \n";
00216 cout << "# Symp. Discr. Alg, p.472 (2002), S.Fortune's Voronoi algorithm and code " ;
00217 #ifndef DROP_CGAL
00218 cout << endl << "# and CGAL: http://www.cgal.org/";
00219 #endif // DROP_CGAL
00220 cout << ".\n";
00221 cout << "#-------------------------------------------------------------------------\n";
00222
00223 cout.flush();
00224 }
00225
00226
00227
00228 void ClusterSequence::_decant_options(const JetDefinition & jet_def,
00229 const bool & writeout_combinations) {
00230
00231
00232 _print_banner();
00233
00234
00235 _jet_def = jet_def;
00236
00237 _writeout_combinations = writeout_combinations;
00238 _jet_algorithm = jet_def.jet_algorithm();
00239 _Rparam = jet_def.R(); _R2 = _Rparam*_Rparam; _invR2 = 1.0/_R2;
00240 _strategy = jet_def.strategy();
00241
00242
00243 _plugin_activated = false;
00244
00245 }
00246
00247
00248
00249
00250 void ClusterSequence::_fill_initial_history () {
00251
00252
00253
00254
00255 _jets.reserve(_jets.size()*2);
00256 _history.reserve(_jets.size()*2);
00257
00258 _Qtot = 0;
00259
00260 for (int i = 0; i < static_cast<int>(_jets.size()) ; i++) {
00261 history_element element;
00262 element.parent1 = InexistentParent;
00263 element.parent2 = InexistentParent;
00264 element.child = Invalid;
00265 element.jetp_index = i;
00266 element.dij = 0.0;
00267 element.max_dij_so_far = 0.0;
00268
00269 _history.push_back(element);
00270
00271
00272 _jet_def.recombiner()->preprocess(_jets[i]);
00273
00274
00275 _jets[i].set_cluster_hist_index(i);
00276
00277
00278 _Qtot += _jets[i].E();
00279 }
00280 _initial_n = _jets.size();
00281 }
00282
00283
00284
00285
00286
00287 string ClusterSequence::strategy_string () const {
00288 string strategy;
00289 switch(_strategy) {
00290 case NlnN:
00291 strategy = "NlnN"; break;
00292 case NlnN3pi:
00293 strategy = "NlnN3pi"; break;
00294 case NlnN4pi:
00295 strategy = "NlnN4pi"; break;
00296 case N2Plain:
00297 strategy = "N2Plain"; break;
00298 case N2Tiled:
00299 strategy = "N2Tiled"; break;
00300 case N2MinHeapTiled:
00301 strategy = "N2MinHeapTiled"; break;
00302 case N2PoorTiled:
00303 strategy = "N2PoorTiled"; break;
00304 case N3Dumb:
00305 strategy = "N3Dumb"; break;
00306 case NlnNCam4pi:
00307 strategy = "NlnNCam4pi"; break;
00308 case NlnNCam2pi2R:
00309 strategy = "NlnNCam2pi2R"; break;
00310 case NlnNCam:
00311 strategy = "NlnNCam"; break;
00312 case plugin_strategy:
00313 strategy = "plugin strategy"; break;
00314 default:
00315 strategy = "Unrecognized";
00316 }
00317 return strategy;
00318 }
00319
00320
00321 double ClusterSequence::jet_scale_for_algorithm(
00322 const PseudoJet & jet) const {
00323 if (_jet_algorithm == kt_algorithm) {return jet.kt2();}
00324 else if (_jet_algorithm == cambridge_algorithm) {return 1.0;}
00325 else if (_jet_algorithm == antikt_algorithm) {
00326 double kt2=jet.kt2();
00327 return kt2 > 1e-300 ? 1.0/kt2 : 1e300;
00328 } else if (_jet_algorithm == genkt_algorithm) {
00329 double kt2 = jet.kt2();
00330 double p = jet_def().extra_param();
00331 if (p <= 0 && kt2 < 1e-300) kt2 = 1e-300;
00332 return pow(kt2, p);
00333 } else if (_jet_algorithm == cambridge_for_passive_algorithm) {
00334 double kt2 = jet.kt2();
00335 double lim = _jet_def.extra_param();
00336 if (kt2 < lim*lim && kt2 != 0.0) {
00337 return 1.0/kt2;
00338 } else {return 1.0;}
00339 } else {throw Error("Unrecognised jet algorithm");}
00340 }
00341
00342
00343
00347 void ClusterSequence::transfer_from_sequence(ClusterSequence & from_seq) {
00348
00349
00350 _jet_def = from_seq._jet_def ;
00351 _writeout_combinations = from_seq._writeout_combinations ;
00352 _initial_n = from_seq._initial_n ;
00353 _Rparam = from_seq._Rparam ;
00354 _R2 = from_seq._R2 ;
00355 _invR2 = from_seq._invR2 ;
00356 _strategy = from_seq._strategy ;
00357 _jet_algorithm = from_seq._jet_algorithm ;
00358 _plugin_activated = from_seq._plugin_activated ;
00359
00360
00361 _jets = from_seq._jets;
00362 _history = from_seq._history;
00363
00364 _extras = from_seq._extras;
00365
00366 }
00367
00368
00369
00370
00371 void ClusterSequence::plugin_record_ij_recombination(
00372 int jet_i, int jet_j, double dij,
00373 const PseudoJet & newjet, int & newjet_k) {
00374
00375 plugin_record_ij_recombination(jet_i, jet_j, dij, newjet_k);
00376
00377
00378 int tmp_index = _jets[newjet_k].cluster_hist_index();
00379 _jets[newjet_k] = newjet;
00380 _jets[newjet_k].set_cluster_hist_index(tmp_index);
00381 }
00382
00383
00384
00385
00386 vector<PseudoJet> ClusterSequence::inclusive_jets (const double & ptmin) const{
00387 double dcut = ptmin*ptmin;
00388 int i = _history.size() - 1;
00389 vector<PseudoJet> jets;
00390 if (_jet_algorithm == kt_algorithm) {
00391 while (i >= 0) {
00392
00393
00394
00395 if (_history[i].max_dij_so_far < dcut) {break;}
00396 if (_history[i].parent2 == BeamJet && _history[i].dij >= dcut) {
00397
00398 int parent1 = _history[i].parent1;
00399 jets.push_back(_jets[_history[parent1].jetp_index]);}
00400 i--;
00401 }
00402 } else if (_jet_algorithm == cambridge_algorithm) {
00403 while (i >= 0) {
00404
00405
00406
00407 if (_history[i].parent2 != BeamJet) {break;}
00408 int parent1 = _history[i].parent1;
00409 const PseudoJet & jet = _jets[_history[parent1].jetp_index];
00410 if (jet.perp2() >= dcut) {jets.push_back(jet);}
00411 i--;
00412 }
00413 } else if (_jet_algorithm == plugin_algorithm
00414 || _jet_algorithm == ee_kt_algorithm
00415 || _jet_algorithm == antikt_algorithm
00416 || _jet_algorithm == genkt_algorithm
00417 || _jet_algorithm == ee_genkt_algorithm
00418 || _jet_algorithm == cambridge_for_passive_algorithm) {
00419
00420
00421
00422 while (i >= 0) {
00423 if (_history[i].parent2 == BeamJet) {
00424 int parent1 = _history[i].parent1;
00425 const PseudoJet & jet = _jets[_history[parent1].jetp_index];
00426 if (jet.perp2() >= dcut) {jets.push_back(jet);}
00427 }
00428 i--;
00429 }
00430 } else {throw Error("cs::inclusive_jets(...): Unrecognized jet algorithm");}
00431 return jets;
00432 }
00433
00434
00435
00436
00437
00438 int ClusterSequence::n_exclusive_jets (const double & dcut) const {
00439
00440
00441
00442 int i = _history.size() - 1;
00443 while (i >= 0) {
00444 if (_history[i].max_dij_so_far <= dcut) {break;}
00445 i--;
00446 }
00447 int stop_point = i + 1;
00448
00449
00450 int njets = 2*_initial_n - stop_point;
00451 return njets;
00452 }
00453
00454
00455
00456
00457 vector<PseudoJet> ClusterSequence::exclusive_jets (const double & dcut) const {
00458 int njets = n_exclusive_jets(dcut);
00459 return exclusive_jets(njets);
00460 }
00461
00462
00463
00464
00465 vector<PseudoJet> ClusterSequence::exclusive_jets (const int & njets) const {
00466
00467
00468
00469 assert (njets <= _initial_n);
00470
00471
00472
00473
00474
00475
00476 if (( _jet_def.jet_algorithm() != kt_algorithm) &&
00477 ( _jet_def.jet_algorithm() != cambridge_algorithm) &&
00478 ( _jet_def.jet_algorithm() != ee_kt_algorithm) &&
00479 (((_jet_def.jet_algorithm() != genkt_algorithm) &&
00480 (_jet_def.jet_algorithm() != ee_genkt_algorithm)) ||
00481 (_jet_def.extra_param() <0)) &&
00482 ((_jet_def.jet_algorithm() != plugin_algorithm) ||
00483 (!_jet_def.plugin()->exclusive_sequence_meaningful())) &&
00484 (_n_exclusive_warnings < 5)) {
00485 _n_exclusive_warnings++;
00486 cerr << "FastJet WARNING: dcut and exclusive jets for jet-finders other than kt should be interpreted with care." << endl;
00487 }
00488
00489
00490
00491
00492
00493 int stop_point = 2*_initial_n - njets;
00494
00495
00496
00497 if (2*_initial_n != static_cast<int>(_history.size())) {
00498 ostringstream err;
00499 err << "2*_initial_n != _history.size() -- this endangers internal assumptions!\n";
00500 throw Error(err.str());
00501
00502 }
00503
00504
00505
00506
00507
00508
00509 vector<PseudoJet> jets;
00510 for (unsigned int i = stop_point; i < _history.size(); i++) {
00511 int parent1 = _history[i].parent1;
00512 if (parent1 < stop_point) {
00513 jets.push_back(_jets[_history[parent1].jetp_index]);
00514 }
00515 int parent2 = _history[i].parent2;
00516 if (parent2 < stop_point && parent2 > 0) {
00517 jets.push_back(_jets[_history[parent2].jetp_index]);
00518 }
00519
00520 }
00521
00522
00523 if (static_cast<int>(jets.size()) != njets) {
00524 ostringstream err;
00525 err << "ClusterSequence::exclusive_jets: size of returned vector ("
00526 <<jets.size()<<") does not coincide with requested number of jets ("
00527 <<njets<<")";
00528 throw Error(err.str());
00529 }
00530
00531 return jets;
00532 }
00533
00534
00537 double ClusterSequence::exclusive_dmerge (const int & njets) const {
00538 assert(njets >= 0);
00539 if (njets >= _initial_n) {return 0.0;}
00540 return _history[2*_initial_n-njets-1].dij;
00541 }
00542
00543
00544
00549 double ClusterSequence::exclusive_dmerge_max (const int & njets) const {
00550 assert(njets >= 0);
00551 if (njets >= _initial_n) {return 0.0;}
00552 return _history[2*_initial_n-njets-1].max_dij_so_far;
00553 }
00554
00555
00556
00560 std::vector<PseudoJet> ClusterSequence::exclusive_subjets
00561 (const PseudoJet & jet, const double & dcut) const {
00562
00563 set<const history_element*> subhist;
00564
00565
00566
00567 get_subhist_set(subhist, jet, dcut, 0);
00568
00569
00570 vector<PseudoJet> subjets;
00571 subjets.reserve(subhist.size());
00572 for (set<const history_element*>::iterator elem = subhist.begin();
00573 elem != subhist.end(); elem++) {
00574 subjets.push_back(_jets[(*elem)->jetp_index]);
00575 }
00576 return subjets;
00577 }
00578
00579
00583 int ClusterSequence::n_exclusive_subjets(const PseudoJet & jet,
00584 const double & dcut) const {
00585 set<const history_element*> subhist;
00586
00587
00588 get_subhist_set(subhist, jet, dcut, 0);
00589 return subhist.size();
00590 }
00591
00592
00596 std::vector<PseudoJet> ClusterSequence::exclusive_subjets
00597 (const PseudoJet & jet, int n) const {
00598
00599 set<const history_element*> subhist;
00600
00601
00602
00603 get_subhist_set(subhist, jet, -1.0, n);
00604
00605
00606 vector<PseudoJet> subjets;
00607 subjets.reserve(subhist.size());
00608 for (set<const history_element*>::iterator elem = subhist.begin();
00609 elem != subhist.end(); elem++) {
00610 subjets.push_back(_jets[(*elem)->jetp_index]);
00611 }
00612 return subjets;
00613 }
00614
00615
00616
00621 double ClusterSequence::exclusive_subdmerge(const PseudoJet & jet, int nsub) const {
00622 set<const history_element*> subhist;
00623
00624
00625
00626 get_subhist_set(subhist, jet, -1.0, nsub);
00627
00628 set<const history_element*>::iterator highest = subhist.end();
00629 highest--;
00632 return (*highest)->dij;
00633 }
00634
00635
00636
00642 double ClusterSequence::exclusive_subdmerge_max(const PseudoJet & jet, int nsub) const {
00643
00644 set<const history_element*> subhist;
00645
00646
00647
00648 get_subhist_set(subhist, jet, -1.0, nsub);
00649
00650 set<const history_element*>::iterator highest = subhist.end();
00651 highest--;
00654 return (*highest)->max_dij_so_far;
00655 }
00656
00657
00658
00659
00666 void ClusterSequence::get_subhist_set(set<const history_element*> & subhist,
00667 const PseudoJet & jet,
00668 double dcut, int maxjet) const {
00669 subhist.clear();
00670 subhist.insert(&(_history[jet.cluster_hist_index()]));
00671
00672
00673 int njet = 1;
00674 while (true) {
00675
00676
00677 set<const history_element*>::iterator highest = subhist.end();
00678 assert (highest != subhist.begin());
00679 highest--;
00680 const history_element* elem = *highest;
00681
00682 if (njet == maxjet) break;
00683
00684 if (elem->parent1 < 0) break;
00685
00686 if (elem->max_dij_so_far <= dcut) break;
00687
00688
00689 subhist.erase(highest);
00690 subhist.insert(&(_history[elem->parent1]));
00691 subhist.insert(&(_history[elem->parent2]));
00692 njet++;
00693 }
00694 }
00695
00696
00697
00698 bool ClusterSequence::object_in_jet(const PseudoJet & object,
00699 const PseudoJet & jet) const {
00700
00701
00702
00703 assert(_potentially_valid(object) && _potentially_valid(jet));
00704
00705 const PseudoJet * this_object = &object;
00706 const PseudoJet * childp;
00707 while(true) {
00708 if (this_object->cluster_hist_index() == jet.cluster_hist_index()) {
00709 return true;
00710 } else if (has_child(*this_object, childp)) {this_object = childp;}
00711 else {return false;}
00712 }
00713 }
00714
00715
00721 bool ClusterSequence::has_parents(const PseudoJet & jet, PseudoJet & parent1,
00722 PseudoJet & parent2) const {
00723
00724 const history_element & hist = _history[jet.cluster_hist_index()];
00725
00726
00727
00728 assert ((hist.parent1 >= 0 && hist.parent2 >= 0) ||
00729 (hist.parent1 < 0 && hist.parent2 < 0));
00730
00731 if (hist.parent1 < 0) {
00732 parent1 = PseudoJet(0.0,0.0,0.0,0.0);
00733 parent2 = parent1;
00734 return false;
00735 } else {
00736 parent1 = _jets[_history[hist.parent1].jetp_index];
00737 parent2 = _jets[_history[hist.parent2].jetp_index];
00738
00739 if (parent1.perp2() < parent2.perp2()) swap(parent1,parent2);
00740 return true;
00741 }
00742 }
00743
00744
00747 bool ClusterSequence::has_child(const PseudoJet & jet, PseudoJet & child) const {
00748
00749
00750
00751
00752
00753
00754
00755
00756
00757
00758 const PseudoJet * childp;
00759 bool res = has_child(jet, childp);
00760 if (res) {
00761 child = *childp;
00762 return true;
00763 } else {
00764 child = PseudoJet(0.0,0.0,0.0,0.0);
00765 return false;
00766 }
00767 }
00768
00769 bool ClusterSequence::has_child(const PseudoJet & jet, const PseudoJet * & childp) const {
00770
00771 const history_element & hist = _history[jet.cluster_hist_index()];
00772
00773
00774
00775
00776 if (hist.child >= 0 && _history[hist.child].jetp_index >= 0) {
00777 childp = &(_jets[_history[hist.child].jetp_index]);
00778 return true;
00779 } else {
00780 childp = NULL;
00781 return false;
00782 }
00783 }
00784
00785
00786
00790 bool ClusterSequence::has_partner(const PseudoJet & jet,
00791 PseudoJet & partner) const {
00792
00793 const history_element & hist = _history[jet.cluster_hist_index()];
00794
00795
00796
00797 if (hist.child >= 0 && _history[hist.child].parent2 >= 0) {
00798 const history_element & child_hist = _history[hist.child];
00799 if (child_hist.parent1 == jet.cluster_hist_index()) {
00800
00801
00802 partner = _jets[_history[child_hist.parent2].jetp_index];
00803 } else {
00804
00805 partner = _jets[_history[child_hist.parent1].jetp_index];
00806 }
00807 return true;
00808 } else {
00809 partner = PseudoJet(0.0,0.0,0.0,0.0);
00810 return false;
00811 }
00812 }
00813
00814
00815
00816
00817 vector<PseudoJet> ClusterSequence::constituents (const PseudoJet & jet) const {
00818 vector<PseudoJet> subjets;
00819 add_constituents(jet, subjets);
00820 return subjets;
00821 }
00822
00823
00832 void ClusterSequence::print_jets_for_root(const std::vector<PseudoJet> & jets,
00833 ostream & ostr) const {
00834 for (unsigned i = 0; i < jets.size(); i++) {
00835 ostr << i << " "
00836 << jets[i].px() << " "
00837 << jets[i].py() << " "
00838 << jets[i].pz() << " "
00839 << jets[i].E() << endl;
00840 vector<PseudoJet> cst = constituents(jets[i]);
00841 for (unsigned j = 0; j < cst.size() ; j++) {
00842 ostr << " " << j << " "
00843 << cst[j].rap() << " "
00844 << cst[j].phi() << " "
00845 << cst[j].perp() << endl;
00846 }
00847 ostr << "#END" << endl;
00848 }
00849 }
00850
00851 void ClusterSequence::print_jets_for_root(const std::vector<PseudoJet> & jets,
00852 const std::string & filename,
00853 const std::string & comment ) const {
00854 std::ofstream ostr(filename.c_str());
00855 if (comment != "") ostr << "# " << comment << endl;
00856 print_jets_for_root(jets, ostr);
00857 }
00858
00859
00860
00861
00862
00863
00864
00865
00866
00867
00868
00869
00870
00871
00872
00877 vector<int> ClusterSequence::particle_jet_indices(
00878 const vector<PseudoJet> & jets) const {
00879
00880 vector<int> indices(n_particles());
00881
00882
00883 for (unsigned ipart = 0; ipart < n_particles(); ipart++)
00884 indices[ipart] = -1;
00885
00886
00887
00888 for (unsigned ijet = 0; ijet < jets.size(); ijet++) {
00889
00890 vector<PseudoJet> jet_constituents(constituents(jets[ijet]));
00891
00892 for (unsigned ip = 0; ip < jet_constituents.size(); ip++) {
00893
00894
00895
00896 unsigned iclust = jet_constituents[ip].cluster_hist_index();
00897 unsigned ipart = history()[iclust].jetp_index;
00898 indices[ipart] = ijet;
00899 }
00900 }
00901
00902 return indices;
00903 }
00904
00905
00906
00907
00908 void ClusterSequence::add_constituents (
00909 const PseudoJet & jet, vector<PseudoJet> & subjet_vector) const {
00910
00911 int i = jet.cluster_hist_index();
00912 int parent1 = _history[i].parent1;
00913 int parent2 = _history[i].parent2;
00914
00915 if (parent1 == InexistentParent) {
00916
00917
00918
00919
00920
00921 subjet_vector.push_back(_jets[i]);
00922 return;
00923 }
00924
00925
00926 add_constituents(_jets[_history[parent1].jetp_index], subjet_vector);
00927
00928
00929 if (parent2 != BeamJet) {
00930 add_constituents(_jets[_history[parent2].jetp_index], subjet_vector);
00931 }
00932 }
00933
00934
00935
00936
00937
00938 void ClusterSequence::_add_step_to_history (
00939 const int & step_number, const int & parent1,
00940 const int & parent2, const int & jetp_index,
00941 const double & dij) {
00942
00943 history_element element;
00944 element.parent1 = parent1;
00945 element.parent2 = parent2;
00946 element.jetp_index = jetp_index;
00947 element.child = Invalid;
00948 element.dij = dij;
00949 element.max_dij_so_far = max(dij,_history[_history.size()-1].max_dij_so_far);
00950 _history.push_back(element);
00951
00952 int local_step = _history.size()-1;
00953 assert(local_step == step_number);
00954
00955 assert(parent1 >= 0);
00956 _history[parent1].child = local_step;
00957 if (parent2 >= 0) {_history[parent2].child = local_step;}
00958
00959
00960 if (jetp_index != Invalid) {
00961 assert(jetp_index >= 0);
00962
00963 _jets[jetp_index].set_cluster_hist_index(local_step);
00964 }
00965
00966 if (_writeout_combinations) {
00967 cout << local_step << ": "
00968 << parent1 << " with " << parent2
00969 << "; y = "<< dij<<endl;
00970 }
00971
00972 }
00973
00974
00975
00976
00977
00978
00979
00980
00981
00982 vector<int> ClusterSequence::unique_history_order() const {
00983
00984
00985
00986
00987
00988 valarray<int> lowest_constituent(_history.size());
00989 int hist_n = _history.size();
00990 lowest_constituent = hist_n;
00991 for (int i = 0; i < hist_n; i++) {
00992
00993 lowest_constituent[i] = min(lowest_constituent[i],i);
00994
00995 if (_history[i].child > 0) lowest_constituent[_history[i].child]
00996 = min(lowest_constituent[_history[i].child],lowest_constituent[i]);
00997 }
00998
00999
01000 valarray<bool> extracted(_history.size()); extracted = false;
01001 vector<int> unique_tree;
01002 unique_tree.reserve(_history.size());
01003
01004
01005 for (unsigned i = 0; i < n_particles(); i++) {
01006 if (!extracted[i]) {
01007 unique_tree.push_back(i);
01008 extracted[i] = true;
01009 _extract_tree_children(i, extracted, lowest_constituent, unique_tree);
01010 }
01011 }
01012
01013 return unique_tree;
01014 }
01015
01016
01017
01018 void ClusterSequence::_extract_tree_children(
01019 int position,
01020 valarray<bool> & extracted,
01021 const valarray<int> & lowest_constituent,
01022 vector<int> & unique_tree) const {
01023 if (!extracted[position]) {
01024
01025
01026 _extract_tree_parents(position,extracted,lowest_constituent,unique_tree);
01027 }
01028
01029
01030 int child = _history[position].child;
01031 if (child >= 0) _extract_tree_children(child,extracted,lowest_constituent,unique_tree);
01032 }
01033
01034
01035
01036
01037 vector<PseudoJet> ClusterSequence::unclustered_particles() const {
01038 vector<PseudoJet> unclustered;
01039 for (unsigned i = 0; i < n_particles() ; i++) {
01040 if (_history[i].child == Invalid)
01041 unclustered.push_back(_jets[_history[i].jetp_index]);
01042 }
01043 return unclustered;
01044 }
01045
01046
01047
01048
01049
01050 void ClusterSequence::_extract_tree_parents(
01051 int position,
01052 valarray<bool> & extracted,
01053 const valarray<int> & lowest_constituent,
01054 vector<int> & unique_tree) const {
01055
01056 if (!extracted[position]) {
01057 int parent1 = _history[position].parent1;
01058 int parent2 = _history[position].parent2;
01059
01060
01061 if (parent1 >= 0 && parent2 >= 0) {
01062 if (lowest_constituent[parent1] > lowest_constituent[parent2])
01063 swap(parent1, parent2);
01064 }
01065
01066 if (parent1 >= 0 && !extracted[parent1])
01067 _extract_tree_parents(parent1,extracted,lowest_constituent,unique_tree);
01068 if (parent2 >= 0 && !extracted[parent2])
01069 _extract_tree_parents(parent2,extracted,lowest_constituent,unique_tree);
01070
01071
01072 unique_tree.push_back(position);
01073 extracted[position] = true;
01074 }
01075 }
01076
01077
01078
01082 void ClusterSequence::_do_ij_recombination_step(
01083 const int & jet_i, const int & jet_j,
01084 const double & dij,
01085 int & newjet_k) {
01086
01087
01088 PseudoJet newjet;
01089 _jet_def.recombiner()->recombine(_jets[jet_i], _jets[jet_j], newjet);
01090 _jets.push_back(newjet);
01091
01092
01093
01094
01095 newjet_k = _jets.size()-1;
01096
01097
01098 int newstep_k = _history.size();
01099
01100 _jets[newjet_k].set_cluster_hist_index(newstep_k);
01101
01102
01103 int hist_i = _jets[jet_i].cluster_hist_index();
01104 int hist_j = _jets[jet_j].cluster_hist_index();
01105
01106 _add_step_to_history(newstep_k, min(hist_i, hist_j), max(hist_i,hist_j),
01107 newjet_k, dij);
01108
01109 }
01110
01111
01112
01115 void ClusterSequence::_do_iB_recombination_step(
01116 const int & jet_i, const double & diB) {
01117
01118 int newstep_k = _history.size();
01119
01120
01121 _add_step_to_history(newstep_k,_jets[jet_i].cluster_hist_index(),BeamJet,
01122 Invalid, diB);
01123
01124 }
01125
01126 FASTJET_END_NAMESPACE
01127