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<cmath>
00038 #include<cstdlib>
00039 #include<cassert>
00040 #include<string>
00041
00042 FASTJET_BEGIN_NAMESPACE
00043
00044 using namespace std;
00045
00047 JetAlgorithm ClusterSequence::_default_jet_algorithm = kt_algorithm;
00048
00049
00050
00051
00052 void ClusterSequence::_initialise_and_run (
00053 const double & R,
00054 const Strategy & strategy,
00055 const bool & writeout_combinations) {
00056
00057 JetDefinition jet_def(_default_jet_algorithm, R, strategy);
00058 _initialise_and_run(jet_def, writeout_combinations);
00059 }
00060
00061
00062
00063 void ClusterSequence::_initialise_and_run (
00064 const JetDefinition & jet_def,
00065 const bool & writeout_combinations) {
00066
00067
00068 _decant_options(jet_def, writeout_combinations);
00069
00070
00071
00072 _fill_initial_history();
00073
00074
00075 if (_jet_algorithm == plugin_algorithm) {
00076 _plugin_activated = true;
00077 _jet_def.plugin()->run_clustering( (*this) );
00078 _plugin_activated = false;
00079 return;
00080 }
00081
00082
00083
00084
00085
00086
00087
00088 if (_strategy == Best) {
00089 int N = _jets.size();
00090 if (N > 6200/pow(_Rparam,2.0)
00091 && jet_def.jet_algorithm() == cambridge_algorithm) {
00092 _strategy = NlnNCam;}
00093 else
00094 #ifndef DROP_CGAL
00095 if (N > 16000/pow(_Rparam,1.15)) {
00096 _strategy = NlnN; }
00097 else
00098 #endif // DROP_CGAL
00099 if (N > 450) {
00100 _strategy = N2MinHeapTiled;
00101 }
00102 else if (N > 55*max(0.5,min(1.0,_Rparam))) {
00103 _strategy = N2Tiled;
00104 } else {
00105 _strategy = N2Plain;
00106 }
00107 }
00108
00109
00110
00111 if (_strategy == NlnN || _strategy == NlnN3pi
00112 || _strategy == NlnN4pi ) {
00113 this->_delaunay_cluster();
00114 } else if (_strategy == N3Dumb ) {
00115 this->_really_dumb_cluster();
00116 } else if (_strategy == N2Tiled) {
00117 this->_faster_tiled_N2_cluster();
00118 } else if (_strategy == N2PoorTiled) {
00119 this->_tiled_N2_cluster();
00120 } else if (_strategy == N2Plain) {
00121 this->_simple_N2_cluster();
00122 } else if (_strategy == N2MinHeapTiled) {
00123 this->_minheap_faster_tiled_N2_cluster();
00124 } else if (_strategy == NlnNCam4pi) {
00125 this->_CP2DChan_cluster();
00126 } else if (_strategy == NlnNCam2pi2R) {
00127 this->_CP2DChan_cluster_2pi2R();
00128 } else if (_strategy == NlnNCam) {
00129 this->_CP2DChan_cluster_2piMultD();
00130 } else {
00131 ostringstream err;
00132 err << "Unrecognised value for strategy: "<<_strategy;
00133 throw Error(err.str());
00134
00135 }
00136 }
00137
00138
00139
00140 bool ClusterSequence::_first_time = true;
00141 int ClusterSequence::_n_exclusive_warnings = 0;
00142
00143
00144
00145
00146 string fastjet_version_string() {
00147 return "FastJet version "+string(fastjet_version);
00148 }
00149
00150
00151
00152
00153 void ClusterSequence::_print_banner() {
00154
00155 if (!_first_time) {return;}
00156 _first_time = false;
00157
00158
00159
00160
00161 cout << "#--------------------------------------------------------------------------\n";
00162 cout << "# FastJet release " << fastjet_version << endl;
00163 cout << "# Written by M. Cacciari, G.P. Salam and G. Soyez \n";
00164 cout << "# http://www.lpthe.jussieu.fr/~salam/fastjet \n";
00165 cout << "# \n";
00166 cout << "# Longitudinally invariant Kt, anti-Kt, and inclusive Cambridge/Aachen \n";
00167 cout << "# clustering using fast geometric algorithms, with area measures and optional\n";
00168 cout << "# external jet-finder plugins. \n";
00169 cout << "# Please cite Phys. Lett. B641 (2006) [hep-ph/0512210] if you use this code.\n";
00170 cout << "# \n";
00171 cout << "# This package uses T.Chan's closest pair algorithm, Proc.13th ACM-SIAM \n";
00172 cout << "# Symp. Discr. Alg, p.472 (2002), S.Fortune's Voronoi algorithm and code " ;
00173 #ifndef DROP_CGAL
00174 cout << endl << "# and CGAL: http://www.cgal.org/";
00175 #endif // DROP_CGAL
00176 cout << ".\n";
00177 cout << "#-------------------------------------------------------------------------\n";
00178 }
00179
00180
00181
00182 void ClusterSequence::_decant_options(const JetDefinition & jet_def,
00183 const bool & writeout_combinations) {
00184
00185
00186 _print_banner();
00187
00188
00189 _jet_def = jet_def;
00190
00191 _writeout_combinations = writeout_combinations;
00192 _jet_algorithm = jet_def.jet_algorithm();
00193 _Rparam = jet_def.R(); _R2 = _Rparam*_Rparam; _invR2 = 1.0/_R2;
00194 _strategy = jet_def.strategy();
00195
00196
00197 _plugin_activated = false;
00198
00199 }
00200
00201
00202
00203
00204 void ClusterSequence::_fill_initial_history () {
00205
00206 if (_jets.size() == 0) {throw Error("Cannot run jet-finder on empty event");}
00207
00208
00209 _jets.reserve(_jets.size()*2);
00210 _history.reserve(_jets.size()*2);
00211
00212 for (int i = 0; i < static_cast<int>(_jets.size()) ; i++) {
00213 history_element element;
00214 element.parent1 = InexistentParent;
00215 element.parent2 = InexistentParent;
00216 element.child = Invalid;
00217 element.jetp_index = i;
00218 element.dij = 0.0;
00219 element.max_dij_so_far = 0.0;
00220
00221 _history.push_back(element);
00222
00223
00224 _jet_def.recombiner()->preprocess(_jets[i]);
00225
00226
00227 _jets[i].set_cluster_hist_index(i);
00228 }
00229 _initial_n = _jets.size();
00230 }
00231
00232
00233
00234
00235
00236 string ClusterSequence::strategy_string () const {
00237 string strategy;
00238 switch(_strategy) {
00239 case NlnN:
00240 strategy = "NlnN"; break;
00241 case NlnN3pi:
00242 strategy = "NlnN3pi"; break;
00243 case NlnN4pi:
00244 strategy = "NlnN4pi"; break;
00245 case N2Plain:
00246 strategy = "N2Plain"; break;
00247 case N2Tiled:
00248 strategy = "N2Tiled"; break;
00249 case N2MinHeapTiled:
00250 strategy = "N2MinHeapTiled"; break;
00251 case N2PoorTiled:
00252 strategy = "N2PoorTiled"; break;
00253 case N3Dumb:
00254 strategy = "N3Dumb"; break;
00255 case NlnNCam4pi:
00256 strategy = "NlnNCam4pi"; break;
00257 case NlnNCam2pi2R:
00258 strategy = "NlnNCam2pi2R"; break;
00259 case NlnNCam:
00260 strategy = "NlnNCam"; break;
00261 case plugin_strategy:
00262 strategy = "plugin strategy"; break;
00263 default:
00264 strategy = "Unrecognized";
00265 }
00266 return strategy;
00267 }
00268
00269
00270 double ClusterSequence::jet_scale_for_algorithm(
00271 const PseudoJet & jet) const {
00272 if (_jet_algorithm == kt_algorithm) {return jet.kt2();}
00273 else if (_jet_algorithm == cambridge_algorithm) {return 1.0;}
00274 else if (_jet_algorithm == antikt_algorithm) {
00275 double kt2=jet.kt2();
00276 return kt2 > 1e-300 ? 1.0/kt2 : 1e300;
00277 } else if (_jet_algorithm == cambridge_for_passive_algorithm) {
00278 double kt2 = jet.kt2();
00279 double lim = _jet_def.extra_param();
00280 if (kt2 < lim*lim && kt2 != 0.0) {
00281 return 1.0/kt2;
00282 } else {return 1.0;}
00283 } else {throw Error("Unrecognised jet finder");}
00284 }
00285
00286
00287
00291 void ClusterSequence::transfer_from_sequence(ClusterSequence & from_seq) {
00292
00293
00294 _jet_def = from_seq._jet_def ;
00295 _writeout_combinations = from_seq._writeout_combinations ;
00296 _initial_n = from_seq._initial_n ;
00297 _Rparam = from_seq._Rparam ;
00298 _R2 = from_seq._R2 ;
00299 _invR2 = from_seq._invR2 ;
00300 _strategy = from_seq._strategy ;
00301 _jet_algorithm = from_seq._jet_algorithm ;
00302 _plugin_activated = from_seq._plugin_activated ;
00303
00304
00305 _jets = from_seq._jets;
00306 _history = from_seq._history;
00307
00308 _extras = from_seq._extras;
00309
00310 }
00311
00312
00313
00314
00315 void ClusterSequence::plugin_record_ij_recombination(
00316 int jet_i, int jet_j, double dij,
00317 const PseudoJet & newjet, int & newjet_k) {
00318
00319 plugin_record_ij_recombination(jet_i, jet_j, dij, newjet_k);
00320
00321
00322 int tmp_index = _jets[newjet_k].cluster_hist_index();
00323 _jets[newjet_k] = newjet;
00324 _jets[newjet_k].set_cluster_hist_index(tmp_index);
00325 }
00326
00327
00328
00329
00330 vector<PseudoJet> ClusterSequence::inclusive_jets (const double & ptmin) const{
00331 double dcut = ptmin*ptmin;
00332 int i = _history.size() - 1;
00333 vector<PseudoJet> jets;
00334 if (_jet_algorithm == kt_algorithm) {
00335 while (i >= 0) {
00336
00337
00338
00339 if (_history[i].max_dij_so_far < dcut) {break;}
00340 if (_history[i].parent2 == BeamJet && _history[i].dij >= dcut) {
00341
00342 int parent1 = _history[i].parent1;
00343 jets.push_back(_jets[_history[parent1].jetp_index]);}
00344 i--;
00345 }
00346 } else if (_jet_algorithm == cambridge_algorithm) {
00347 while (i >= 0) {
00348
00349
00350
00351 if (_history[i].parent2 != BeamJet) {break;}
00352 int parent1 = _history[i].parent1;
00353 const PseudoJet & jet = _jets[_history[parent1].jetp_index];
00354 if (jet.perp2() >= dcut) {jets.push_back(jet);}
00355 i--;
00356 }
00357 } else if (_jet_algorithm == plugin_algorithm
00358 || _jet_algorithm == antikt_algorithm
00359 || _jet_algorithm == cambridge_for_passive_algorithm) {
00360
00361
00362
00363 while (i >= 0) {
00364 if (_history[i].parent2 == BeamJet) {
00365 int parent1 = _history[i].parent1;
00366 const PseudoJet & jet = _jets[_history[parent1].jetp_index];
00367 if (jet.perp2() >= dcut) {jets.push_back(jet);}
00368 }
00369 i--;
00370 }
00371 } else {throw Error("Unrecognized jet algorithm");}
00372 return jets;
00373 }
00374
00375
00376
00377
00378
00379 int ClusterSequence::n_exclusive_jets (const double & dcut) const {
00380
00381
00382
00383 int i = _history.size() - 1;
00384 while (i >= 0) {
00385 if (_history[i].max_dij_so_far <= dcut) {break;}
00386 i--;
00387 }
00388 int stop_point = i + 1;
00389
00390
00391 int njets = 2*_initial_n - stop_point;
00392 return njets;
00393 }
00394
00395
00396
00397
00398 vector<PseudoJet> ClusterSequence::exclusive_jets (const double & dcut) const {
00399 int njets = n_exclusive_jets(dcut);
00400 return exclusive_jets(njets);
00401 }
00402
00403
00404
00405
00406 vector<PseudoJet> ClusterSequence::exclusive_jets (const int & njets) const {
00407
00408
00409
00410 assert (njets <= _initial_n);
00411
00412
00413 if (_jet_def.jet_algorithm() != kt_algorithm && _n_exclusive_warnings < 5) {
00414 _n_exclusive_warnings++;
00415 cerr << "FastJet WARNING: dcut and exclusive jets for jet-finders other than kt should be interpreted with care." << endl;
00416 }
00417
00418
00419
00420
00421
00422 int stop_point = 2*_initial_n - njets;
00423
00424
00425
00426 if (2*_initial_n != static_cast<int>(_history.size())) {
00427 ostringstream err;
00428 err << "2*_initial_n != _history.size() -- this endangers internal assumptions!\n";
00429 throw Error(err.str());
00430
00431 }
00432
00433
00434
00435
00436
00437
00438 vector<PseudoJet> jets;
00439 for (unsigned int i = stop_point; i < _history.size(); i++) {
00440 int parent1 = _history[i].parent1;
00441 if (parent1 < stop_point) {
00442 jets.push_back(_jets[_history[parent1].jetp_index]);
00443 }
00444 int parent2 = _history[i].parent2;
00445 if (parent2 < stop_point && parent2 > 0) {
00446 jets.push_back(_jets[_history[parent2].jetp_index]);
00447 }
00448
00449 }
00450
00451
00452 if (static_cast<int>(jets.size()) != njets) {
00453 ostringstream err;
00454 err << "ClusterSequence::exclusive_jets: size of returned vector ("
00455 <<jets.size()<<") does not coincide with requested number of jets ("
00456 <<njets<<")";
00457 throw Error(err.str());
00458 }
00459
00460 return jets;
00461 }
00462
00463
00466 double ClusterSequence::exclusive_dmerge (const int & njets) const {
00467 assert(njets > 0);
00468 if (njets >= _initial_n) {return 0.0;}
00469 return _history[2*_initial_n-njets-1].dij;
00470 }
00471
00472
00473
00478 double ClusterSequence::exclusive_dmerge_max (const int & njets) const {
00479 assert(njets > 0);
00480 if (njets >= _initial_n) {return 0.0;}
00481 return _history[2*_initial_n-njets-1].max_dij_so_far;
00482 }
00483
00484
00485
00486
00487 bool ClusterSequence::object_in_jet(const PseudoJet & object,
00488 const PseudoJet & jet) const {
00489
00490
00491
00492 assert(_potentially_valid(object) && _potentially_valid(jet));
00493
00494 const PseudoJet * this_object = &object;
00495 const PseudoJet * childp;
00496 while(true) {
00497 if (this_object->cluster_hist_index() == jet.cluster_hist_index()) {
00498 return true;
00499 } else if (has_child(*this_object, childp)) {this_object = childp;}
00500 else {return false;}
00501 }
00502 }
00503
00504
00510 bool ClusterSequence::has_parents(const PseudoJet & jet, PseudoJet & parent1,
00511 PseudoJet & parent2) const {
00512
00513 const history_element & hist = _history[jet.cluster_hist_index()];
00514
00515
00516
00517 assert ((hist.parent1 >= 0 && hist.parent2 >= 0) ||
00518 (hist.parent1 < 0 && hist.parent2 < 0));
00519
00520 if (hist.parent1 < 0) {
00521 parent1 = PseudoJet(0.0,0.0,0.0,0.0);
00522 parent2 = parent1;
00523 return false;
00524 } else {
00525 parent1 = _jets[_history[hist.parent1].jetp_index];
00526 parent2 = _jets[_history[hist.parent2].jetp_index];
00527
00528 if (parent1.perp2() < parent2.perp2()) swap(parent1,parent2);
00529 return true;
00530 }
00531 }
00532
00533
00536 bool ClusterSequence::has_child(const PseudoJet & jet, PseudoJet & child) const {
00537
00538
00539
00540
00541
00542
00543
00544
00545
00546
00547 const PseudoJet * childp;
00548 bool res = has_child(jet, childp);
00549 if (res) {
00550 child = *childp;
00551 return true;
00552 } else {
00553 child = PseudoJet(0.0,0.0,0.0,0.0);
00554 return false;
00555 }
00556 }
00557
00558 bool ClusterSequence::has_child(const PseudoJet & jet, const PseudoJet * & childp) const {
00559
00560 const history_element & hist = _history[jet.cluster_hist_index()];
00561
00562
00563
00564
00565 if (hist.child >= 0 && _history[hist.child].jetp_index >= 0) {
00566 childp = &(_jets[_history[hist.child].jetp_index]);
00567 return true;
00568 } else {
00569 childp = NULL;
00570 return false;
00571 }
00572 }
00573
00574
00575
00579 bool ClusterSequence::has_partner(const PseudoJet & jet,
00580 PseudoJet & partner) const {
00581
00582 const history_element & hist = _history[jet.cluster_hist_index()];
00583
00584
00585
00586 if (hist.child >= 0 && _history[hist.child].parent2 >= 0) {
00587 const history_element & child_hist = _history[hist.child];
00588 if (child_hist.parent1 == jet.cluster_hist_index()) {
00589
00590
00591 partner = _jets[_history[child_hist.parent2].jetp_index];
00592 } else {
00593
00594 partner = _jets[_history[child_hist.parent1].jetp_index];
00595 }
00596 return true;
00597 } else {
00598 partner = PseudoJet(0.0,0.0,0.0,0.0);
00599 return false;
00600 }
00601 }
00602
00603
00604
00605
00606 vector<PseudoJet> ClusterSequence::constituents (const PseudoJet & jet) const {
00607 vector<PseudoJet> subjets;
00608 add_constituents(jet, subjets);
00609 return subjets;
00610 }
00611
00612
00621 void ClusterSequence::print_jets_for_root(const std::vector<PseudoJet> & jets,
00622 ostream & ostr) const {
00623 for (unsigned i = 0; i < jets.size(); i++) {
00624 ostr << i << " "
00625 << jets[i].px() << " "
00626 << jets[i].py() << " "
00627 << jets[i].pz() << " "
00628 << jets[i].E() << endl;
00629 vector<PseudoJet> cst = constituents(jets[i]);
00630 for (unsigned j = 0; j < cst.size() ; j++) {
00631 ostr << " " << j << " "
00632 << cst[j].rap() << " "
00633 << cst[j].phi() << " "
00634 << cst[j].perp() << endl;
00635 }
00636 ostr << "#END" << endl;
00637 }
00638 }
00639
00640
00641
00642
00643
00644
00645
00646
00647
00648
00649
00650
00651
00652
00657 vector<int> ClusterSequence::particle_jet_indices(
00658 const vector<PseudoJet> & jets) const {
00659
00660 vector<int> indices(n_particles());
00661
00662
00663 for (unsigned ipart = 0; ipart < n_particles(); ipart++)
00664 indices[ipart] = -1;
00665
00666
00667
00668 for (unsigned ijet = 0; ijet < jets.size(); ijet++) {
00669
00670 vector<PseudoJet> jet_constituents(constituents(jets[ijet]));
00671
00672 for (unsigned ip = 0; ip < jet_constituents.size(); ip++) {
00673
00674
00675
00676 unsigned iclust = jet_constituents[ip].cluster_hist_index();
00677 unsigned ipart = history()[iclust].jetp_index;
00678 indices[ipart] = ijet;
00679 }
00680 }
00681
00682 return indices;
00683 }
00684
00685
00686
00687
00688 void ClusterSequence::add_constituents (
00689 const PseudoJet & jet, vector<PseudoJet> & subjet_vector) const {
00690
00691 int i = jet.cluster_hist_index();
00692 int parent1 = _history[i].parent1;
00693 int parent2 = _history[i].parent2;
00694
00695 if (parent1 == InexistentParent) {
00696
00697
00698 subjet_vector.push_back(jet);
00699 return;
00700 }
00701
00702
00703 add_constituents(_jets[_history[parent1].jetp_index], subjet_vector);
00704
00705
00706 if (parent2 != BeamJet) {
00707 add_constituents(_jets[_history[parent2].jetp_index], subjet_vector);
00708 }
00709 }
00710
00711
00712
00713
00714
00715 void ClusterSequence::_add_step_to_history (
00716 const int & step_number, const int & parent1,
00717 const int & parent2, const int & jetp_index,
00718 const double & dij) {
00719
00720 history_element element;
00721 element.parent1 = parent1;
00722 element.parent2 = parent2;
00723 element.jetp_index = jetp_index;
00724 element.child = Invalid;
00725 element.dij = dij;
00726 element.max_dij_so_far = max(dij,_history[_history.size()-1].max_dij_so_far);
00727 _history.push_back(element);
00728
00729 int local_step = _history.size()-1;
00730 assert(local_step == step_number);
00731
00732 assert(parent1 >= 0);
00733 _history[parent1].child = local_step;
00734 if (parent2 >= 0) {_history[parent2].child = local_step;}
00735
00736
00737 if (jetp_index != Invalid) {
00738 assert(jetp_index >= 0);
00739
00740 _jets[jetp_index].set_cluster_hist_index(local_step);
00741 }
00742
00743 if (_writeout_combinations) {
00744 cout << local_step << ": "
00745 << parent1 << " with " << parent2
00746 << "; y = "<< dij<<endl;
00747 }
00748
00749 }
00750
00751
00752
00753
00754
00755
00756
00757
00758
00759 vector<int> ClusterSequence::unique_history_order() const {
00760
00761
00762
00763
00764
00765 valarray<int> lowest_constituent(_history.size());
00766 int hist_n = _history.size();
00767 lowest_constituent = hist_n;
00768 for (int i = 0; i < hist_n; i++) {
00769
00770 lowest_constituent[i] = min(lowest_constituent[i],i);
00771
00772 if (_history[i].child > 0) lowest_constituent[_history[i].child]
00773 = min(lowest_constituent[_history[i].child],lowest_constituent[i]);
00774 }
00775
00776
00777 valarray<bool> extracted(_history.size()); extracted = false;
00778 vector<int> unique_tree;
00779 unique_tree.reserve(_history.size());
00780
00781
00782 for (unsigned i = 0; i < n_particles(); i++) {
00783 if (!extracted[i]) {
00784 unique_tree.push_back(i);
00785 extracted[i] = true;
00786 _extract_tree_children(i, extracted, lowest_constituent, unique_tree);
00787 }
00788 }
00789
00790 return unique_tree;
00791 }
00792
00793
00794
00795 void ClusterSequence::_extract_tree_children(
00796 int position,
00797 valarray<bool> & extracted,
00798 const valarray<int> & lowest_constituent,
00799 vector<int> & unique_tree) const {
00800 if (!extracted[position]) {
00801
00802
00803 _extract_tree_parents(position,extracted,lowest_constituent,unique_tree);
00804 }
00805
00806
00807 int child = _history[position].child;
00808 if (child >= 0) _extract_tree_children(child,extracted,lowest_constituent,unique_tree);
00809 }
00810
00811
00812
00813
00814 vector<PseudoJet> ClusterSequence::unclustered_particles() const {
00815 vector<PseudoJet> unclustered;
00816 for (unsigned i = 0; i < n_particles() ; i++) {
00817 if (_history[i].child == Invalid)
00818 unclustered.push_back(_jets[_history[i].jetp_index]);
00819 }
00820 return unclustered;
00821 }
00822
00823
00824
00825
00826
00827 void ClusterSequence::_extract_tree_parents(
00828 int position,
00829 valarray<bool> & extracted,
00830 const valarray<int> & lowest_constituent,
00831 vector<int> & unique_tree) const {
00832
00833 if (!extracted[position]) {
00834 int parent1 = _history[position].parent1;
00835 int parent2 = _history[position].parent2;
00836
00837
00838 if (parent1 >= 0 && parent2 >= 0) {
00839 if (lowest_constituent[parent1] > lowest_constituent[parent2])
00840 swap(parent1, parent2);
00841 }
00842
00843 if (parent1 >= 0 && !extracted[parent1])
00844 _extract_tree_parents(parent1,extracted,lowest_constituent,unique_tree);
00845 if (parent2 >= 0 && !extracted[parent2])
00846 _extract_tree_parents(parent2,extracted,lowest_constituent,unique_tree);
00847
00848
00849 unique_tree.push_back(position);
00850 extracted[position] = true;
00851 }
00852 }
00853
00854
00855
00859 void ClusterSequence::_do_ij_recombination_step(
00860 const int & jet_i, const int & jet_j,
00861 const double & dij,
00862 int & newjet_k) {
00863
00864
00865 PseudoJet newjet;
00866 _jet_def.recombiner()->recombine(_jets[jet_i], _jets[jet_j], newjet);
00867 _jets.push_back(newjet);
00868
00869
00870
00871
00872 newjet_k = _jets.size()-1;
00873
00874
00875 int newstep_k = _history.size();
00876
00877 _jets[newjet_k].set_cluster_hist_index(newstep_k);
00878
00879
00880 int hist_i = _jets[jet_i].cluster_hist_index();
00881 int hist_j = _jets[jet_j].cluster_hist_index();
00882
00883 _add_step_to_history(newstep_k, min(hist_i, hist_j), max(hist_i,hist_j),
00884 newjet_k, dij);
00885
00886 }
00887
00888
00889
00892 void ClusterSequence::_do_iB_recombination_step(
00893 const int & jet_i, const double & diB) {
00894
00895 int newstep_k = _history.size();
00896
00897
00898 _add_step_to_history(newstep_k,_jets[jet_i].cluster_hist_index(),BeamJet,
00899 Invalid, diB);
00900
00901 }
00902
00903 FASTJET_END_NAMESPACE
00904