31 #include "fastjet/Error.hh"
32 #include "fastjet/PseudoJet.hh"
33 #include "fastjet/ClusterSequence.hh"
34 #include "fastjet/ClusterSequenceStructure.hh"
35 #include "fastjet/version.hh"
36 #include "fastjet/internal/LazyTiling9Alt.hh"
37 #include "fastjet/internal/LazyTiling9.hh"
38 #include "fastjet/internal/LazyTiling25.hh"
40 #include "fastjet/internal/LazyTiling9SeparateGhosts.hh"
51 FASTJET_BEGIN_NAMESPACE
147 #ifdef FASTJET_HAVE_LIMITED_THREAD_SAFETY
148 atomic<ostream *> ClusterSequence::_fastjet_banner_ostr{&cout};
150 ostream * ClusterSequence::_fastjet_banner_ostr = &cout;
151 #endif // FASTJET_HAVE_LIMITED_THREAD_SAFETY
155 ClusterSequence::~ClusterSequence () {
158 if (_structure_shared_ptr){
159 ClusterSequenceStructure* csi =
dynamic_cast<ClusterSequenceStructure*
>(_structure_shared_ptr.get());
165 csi->set_associated_cs(NULL);
177 if (_deletes_self_when_unused) {
178 _structure_shared_ptr.set_count(_structure_shared_ptr.use_count()
179 + _structure_use_count_after_construction);
216 void ClusterSequence::signal_imminent_self_deletion()
const {
230 assert(_deletes_self_when_unused);
231 _deletes_self_when_unused =
false;
248 void ClusterSequence::_initialise_and_run (
250 const bool & writeout_combinations) {
253 _decant_options(jet_def_in, writeout_combinations);
256 _initialise_and_run_no_decant();
260 void ClusterSequence::_initialise_and_run_no_decant () {
264 _fill_initial_history();
267 if (n_particles() == 0)
return;
270 if (_jet_algorithm == plugin_algorithm) {
272 _plugin_activated =
true;
274 _jet_def.plugin()->run_clustering( (*
this) );
275 _plugin_activated =
false;
276 _update_structure_use_count();
278 }
else if (_jet_algorithm == ee_kt_algorithm ||
279 _jet_algorithm == ee_genkt_algorithm) {
282 if (_jet_algorithm == ee_kt_algorithm) {
286 assert(_Rparam > 2.0);
300 _R2 = 2 * ( 3.0 + cos(_Rparam) );
302 _R2 = 2 * ( 1.0 - cos(_Rparam) );
306 _simple_N2_cluster_EEBriefJet();
308 }
else if (_jet_algorithm == undefined_jet_algorithm) {
309 throw Error(
"A ClusterSequence cannot be created with an uninitialised JetDefinition");
323 if (_strategy == Best) {
324 _strategy = _best_strategy();
329 }
else if (_strategy == BestFJ30) {
330 int N = _jets.size();
336 if (min(1.0,max(0.1,_Rparam)*3.3)*N <= 30) {
338 }
else if (N > 6200/pow(_Rparam,2.0) && _jet_def.jet_algorithm() == cambridge_algorithm) {
341 }
else if ((N > 16000/pow(_Rparam,1.15) && _jet_def.jet_algorithm() != antikt_algorithm)
342 || N > 35000/pow(_Rparam,1.15)) {
345 }
else if (N <= 450) {
356 if (_Rparam >= twopi) {
357 if ( _strategy == NlnN
358 || _strategy == NlnN3pi
359 || _strategy == NlnNCam
360 || _strategy == NlnNCam2pi2R
361 || _strategy == NlnNCam4pi) {
368 if (_jet_def.strategy() != Best && _strategy != _jet_def.strategy()) {
370 oss <<
"Cluster strategy " << strategy_string(_jet_def.strategy())
371 <<
" automatically changed to " << strategy_string()
372 <<
" because the former is not supported for R = " << _Rparam
374 _changed_strategy_warning.warn(oss.str());
384 if (_strategy == N2Plain) {
386 this->_simple_N2_cluster_BriefJet();
387 }
else if (_strategy == N2Tiled) {
388 this->_faster_tiled_N2_cluster();
389 }
else if (_strategy == N2MinHeapTiled) {
390 this->_minheap_faster_tiled_N2_cluster();
391 }
else if (_strategy == N2MHTLazy9Alt) {
394 _plugin_activated =
true;
395 LazyTiling9Alt tiling(*
this);
397 _plugin_activated =
false;
399 }
else if (_strategy == N2MHTLazy25) {
402 _plugin_activated =
true;
403 LazyTiling25 tiling(*
this);
405 _plugin_activated =
false;
407 }
else if (_strategy == N2MHTLazy9) {
410 _plugin_activated =
true;
411 LazyTiling9 tiling(*
this);
413 _plugin_activated =
false;
415 }
else if (_strategy == N2MHTLazy9AntiKtSeparateGhosts) {
419 _plugin_activated =
true;
420 LazyTiling9SeparateGhosts tiling(*
this);
422 _plugin_activated =
false;
424 throw Error(
"N2MHTLazy9AntiKtSeparateGhosts strategy not supported with FJCORE");
427 }
else if (_strategy == NlnN) {
428 this->_delaunay_cluster();
429 }
else if (_strategy == NlnNCam) {
430 this->_CP2DChan_cluster_2piMultD();
431 }
else if (_strategy == NlnN3pi || _strategy == NlnN4pi ) {
432 this->_delaunay_cluster();
433 }
else if (_strategy == N3Dumb ) {
434 this->_really_dumb_cluster();
435 }
else if (_strategy == N2PoorTiled) {
436 this->_tiled_N2_cluster();
437 }
else if (_strategy == NlnNCam4pi) {
438 this->_CP2DChan_cluster();
439 }
else if (_strategy == NlnNCam2pi2R) {
440 this->_CP2DChan_cluster_2pi2R();
443 err <<
"Unrecognised value for strategy: "<<_strategy;
444 throw Error(err.str());
451 thread_safety_helpers::FirstTimeTrue ClusterSequence::_first_time;
452 LimitedWarning ClusterSequence::_exclusive_warnings;
458 return "FastJet version "+string(fastjet_version);
464 void ClusterSequence::print_banner() {
466 if (!_first_time())
return;
469 ostream * ostr = _fastjet_banner_ostr;
472 (*ostr) <<
"#--------------------------------------------------------------------------\n";
473 (*ostr) <<
"# FastJet release " << fastjet_version << endl;
474 (*ostr) <<
"# M. Cacciari, G.P. Salam and G. Soyez \n";
475 (*ostr) <<
"# A software package for jet finding and analysis at colliders \n";
476 (*ostr) <<
"# http://fastjet.fr \n";
478 (*ostr) <<
"# Please cite EPJC72(2012)1896 [arXiv:1111.6097] if you use this package\n";
479 (*ostr) <<
"# for scientific work and optionally PLB641(2006)57 [hep-ph/0512210]. \n";
481 (*ostr) <<
"# FastJet is provided without warranty under the GNU GPL v2 or higher. \n";
482 (*ostr) <<
"# It uses T. Chan's closest pair algorithm, S. Fortune's Voronoi code";
484 (*ostr) <<
",\n# CGAL ";
488 (*ostr) <<
"and 3rd party plugin jet algorithms. See COPYING file for details.\n";
489 (*ostr) <<
"#--------------------------------------------------------------------------\n";
497 const bool & writeout_combinations) {
499 _jet_def = jet_def_in;
500 _writeout_combinations = writeout_combinations;
504 _decant_options_partial();
509 void ClusterSequence::_decant_options_partial() {
513 _jet_algorithm = _jet_def.jet_algorithm();
514 _Rparam = _jet_def.R(); _R2 = _Rparam*_Rparam; _invR2 = 1.0/_R2;
515 _strategy = _jet_def.strategy();
518 _plugin_activated =
false;
522 _update_structure_use_count();
528 void ClusterSequence::_fill_initial_history () {
533 _jets.reserve(_jets.size()*2);
534 _history.reserve(_jets.size()*2);
538 for (
int i = 0; i < static_cast<int>(_jets.size()) ; i++) {
540 element.parent1 = InexistentParent;
541 element.
parent2 = InexistentParent;
542 element.
child = Invalid;
547 _history.push_back(element);
550 _jet_def.recombiner()->preprocess(_jets[i]);
553 _jets[i].set_cluster_hist_index(i);
554 _set_structure_shared_ptr(_jets[i]);
557 _Qtot += _jets[i].E();
559 _initial_n = _jets.size();
560 _deletes_self_when_unused =
false;
565 string ClusterSequence::strategy_string (
Strategy strategy_in)
const {
567 switch(strategy_in) {
569 strategy =
"NlnN";
break;
571 strategy =
"NlnN3pi";
break;
573 strategy =
"NlnN4pi";
break;
575 strategy =
"N2Plain";
break;
577 strategy =
"N2Tiled";
break;
579 strategy =
"N2MinHeapTiled";
break;
581 strategy =
"N2PoorTiled";
break;
583 strategy =
"N2MHTLazy9";
break;
585 strategy =
"N2MHTLazy9Alt";
break;
587 strategy =
"N2MHTLazy25";
break;
589 strategy =
"N2MHTLazy9AntiKtSeparateGhosts";
break;
591 strategy =
"N3Dumb";
break;
593 strategy =
"NlnNCam4pi";
break;
595 strategy =
"NlnNCam2pi2R";
break;
597 strategy =
"NlnNCam";
break;
599 strategy =
"plugin strategy";
break;
601 strategy =
"Unrecognized";
607 double ClusterSequence::jet_scale_for_algorithm(
612 double kt2=jet.
kt2();
613 return kt2 > 1e-300 ? 1.0/kt2 : 1e300;
615 double kt2 = jet.
kt2();
616 double p = jet_def().extra_param();
617 if (p <= 0 && kt2 < 1e-300) kt2 = 1e-300;
620 double kt2 = jet.
kt2();
621 double lim = _jet_def.extra_param();
622 if (kt2 < lim*lim && kt2 != 0.0) {
625 }
else {
throw Error(
"Unrecognised jet algorithm");}
644 int N = _jets.size();
647 double bounded_R = max(_Rparam, 0.1);
651 if (N <= 30 || N <= 39.0/(bounded_R + 0.6)) {
662 const static _Parabola N_Tiled_to_MHT_lowR (-45.4947,54.3528,44.6283);
663 const static _Parabola L_MHT_to_MHTLazy9_lowR (0.677807,-1.05006,10.6994);
664 const static _Parabola L_MHTLazy9_to_MHTLazy25_akt_lowR(0.169967,-0.512589,12.1572);
665 const static _Parabola L_MHTLazy9_to_MHTLazy25_kt_lowR (0.16237,-0.484612,12.3373);
666 const static _Parabola L_MHTLazy9_to_MHTLazy25_cam_lowR = L_MHTLazy9_to_MHTLazy25_kt_lowR;
667 const static _Parabola L_MHTLazy25_to_NlnN_akt_lowR (0.0472051,-0.22043,15.9196);
668 const static _Parabola L_MHTLazy25_to_NlnN_kt_lowR (0.118609,-0.326811,14.8287);
669 const static _Parabola L_MHTLazy25_to_NlnN_cam_lowR (0.10119,-0.295748,14.3924);
671 const static _Line L_Tiled_to_MHTLazy9_medR (-1.31304,7.29621);
672 const static _Parabola L_MHTLazy9_to_MHTLazy25_akt_medR = L_MHTLazy9_to_MHTLazy25_akt_lowR;
673 const static _Parabola L_MHTLazy9_to_MHTLazy25_kt_medR = L_MHTLazy9_to_MHTLazy25_kt_lowR;
674 const static _Parabola L_MHTLazy9_to_MHTLazy25_cam_medR = L_MHTLazy9_to_MHTLazy25_cam_lowR;
675 const static _Parabola L_MHTLazy25_to_NlnN_akt_medR = L_MHTLazy25_to_NlnN_akt_lowR;
676 const static _Parabola L_MHTLazy25_to_NlnN_kt_medR = L_MHTLazy25_to_NlnN_kt_lowR;
677 const static _Parabola L_MHTLazy25_to_NlnN_cam_medR = L_MHTLazy25_to_NlnN_cam_lowR;
679 const static double N_Plain_to_MHTLazy9_largeR = 75;
680 const static double N_MHTLazy9_to_MHTLazy25_akt_largeR = 700;
681 const static double N_MHTLazy9_to_MHTLazy25_kt_largeR = 1000;
682 const static double N_MHTLazy9_to_MHTLazy25_cam_largeR = 1000;
683 const static double N_MHTLazy25_to_NlnN_akt_largeR = 100000;
684 const static double N_MHTLazy25_to_NlnN_kt_largeR = 40000;
685 const static double N_MHTLazy25_to_NlnN_cam_largeR = 15000;
696 double p = jet_def().extra_param();
704 jet_algorithm = _jet_algorithm;
707 if (bounded_R < 0.65) {
709 if (N < N_Tiled_to_MHT_lowR(bounded_R))
return N2Tiled;
710 double logN = log(
double(N));
711 if (logN < L_MHT_to_MHTLazy9_lowR(bounded_R))
return N2MinHeapTiled;
714 if (logN < L_MHTLazy9_to_MHTLazy25_akt_lowR(bounded_R))
return N2MHTLazy9;
715 else if (logN < L_MHTLazy25_to_NlnN_akt_lowR(bounded_R))
return N2MHTLazy25;
718 if (logN < L_MHTLazy9_to_MHTLazy25_kt_lowR(bounded_R))
return N2MHTLazy9;
719 else if (logN < L_MHTLazy25_to_NlnN_kt_lowR(bounded_R))
return N2MHTLazy25;
722 if (logN < L_MHTLazy9_to_MHTLazy25_cam_lowR(bounded_R))
return N2MHTLazy9;
723 else if (logN < L_MHTLazy25_to_NlnN_cam_lowR(bounded_R))
return N2MHTLazy25;
727 }
else if (bounded_R < 0.5*pi) {
729 double logN = log(
double(N));
730 if (logN < L_Tiled_to_MHTLazy9_medR(bounded_R))
return N2Tiled;
733 if (logN < L_MHTLazy9_to_MHTLazy25_akt_medR(bounded_R))
return N2MHTLazy9;
734 else if (logN < L_MHTLazy25_to_NlnN_akt_medR(bounded_R))
return N2MHTLazy25;
737 if (logN < L_MHTLazy9_to_MHTLazy25_kt_medR(bounded_R))
return N2MHTLazy9;
738 else if (logN < L_MHTLazy25_to_NlnN_kt_medR(bounded_R))
return N2MHTLazy25;
741 if (logN < L_MHTLazy9_to_MHTLazy25_cam_medR(bounded_R))
return N2MHTLazy9;
742 else if (logN < L_MHTLazy25_to_NlnN_cam_medR(bounded_R))
return N2MHTLazy25;
748 if (N < N_Plain_to_MHTLazy9_largeR)
return N2Plain;
751 if (N < N_MHTLazy9_to_MHTLazy25_akt_largeR)
return N2MHTLazy9;
752 else if (N < N_MHTLazy25_to_NlnN_akt_largeR)
return N2MHTLazy25;
755 if (N < N_MHTLazy9_to_MHTLazy25_kt_largeR)
return N2MHTLazy9;
756 else if (N < N_MHTLazy25_to_NlnN_kt_largeR)
return N2MHTLazy25;
759 if (N < N_MHTLazy9_to_MHTLazy25_cam_largeR)
return N2MHTLazy9;
760 else if (N < N_MHTLazy25_to_NlnN_cam_largeR)
return N2MHTLazy25;
769 assert(0 &&
"Code should never reach here");
823 _deletes_self_when_unused =
false;
824 transfer_from_sequence(cs);
842 if (will_delete_self_when_unused())
843 throw(
Error(
"cannot use CS::transfer_from_sequence after a call to delete_self_when_unused()"));
846 _jet_def = from_seq._jet_def ;
847 _writeout_combinations = from_seq._writeout_combinations ;
848 _initial_n = from_seq._initial_n ;
849 _Rparam = from_seq._Rparam ;
851 _invR2 = from_seq._invR2 ;
852 _strategy = from_seq._strategy ;
853 _jet_algorithm = from_seq._jet_algorithm ;
854 _plugin_activated = from_seq._plugin_activated ;
860 _jets = (*action_on_jets)(from_seq.
_jets);
862 _jets = from_seq.
_jets;
866 _extras = from_seq._extras;
869 if (_structure_shared_ptr) {
873 if (_deletes_self_when_unused)
throw Error(
"transfer_from_sequence cannot be used for a cluster sequence that deletes self when unused");
884 _update_structure_use_count();
886 for (
unsigned int i=0; i<_jets.size(); i++){
889 _jets[i].set_cluster_hist_index(from_seq.
_jets[i].cluster_hist_index());
892 _set_structure_shared_ptr(_jets[i]);
900 void ClusterSequence::plugin_record_ij_recombination(
901 int jet_i,
int jet_j,
double dij,
902 const PseudoJet & newjet,
int & newjet_k) {
904 plugin_record_ij_recombination(jet_i, jet_j, dij, newjet_k);
907 int tmp_index = _jets[newjet_k].cluster_hist_index();
908 _jets[newjet_k] = newjet;
910 _set_structure_shared_ptr(_jets[newjet_k]);
916 vector<PseudoJet> ClusterSequence::inclusive_jets (
const double ptmin)
const{
917 double dcut = ptmin*ptmin;
918 int i = _history.size() - 1;
919 vector<PseudoJet> jets_local;
925 if (_history[i].max_dij_so_far < dcut) {
break;}
926 if (_history[i].parent2 == BeamJet && _history[i].dij >= dcut) {
928 int parent1 = _history[i].parent1;
929 jets_local.push_back(_jets[_history[parent1].jetp_index]);}
937 if (_history[i].parent2 != BeamJet) {
break;}
938 int parent1 = _history[i].parent1;
939 const PseudoJet & jet = _jets[_history[parent1].jetp_index];
940 if (jet.
perp2() >= dcut) {jets_local.push_back(jet);}
953 if (_history[i].parent2 == BeamJet) {
954 int parent1 = _history[i].parent1;
955 const PseudoJet & jet = _jets[_history[parent1].jetp_index];
956 if (jet.
perp2() >= dcut) {jets_local.push_back(jet);}
960 }
else {
throw Error(
"cs::inclusive_jets(...): Unrecognized jet algorithm");}
968 int ClusterSequence::n_exclusive_jets (
const double dcut)
const {
972 int i = _history.size() - 1;
974 if (_history[i].max_dij_so_far <= dcut) {
break;}
977 int stop_point = i + 1;
980 int njets = 2*_initial_n - stop_point;
987 vector<PseudoJet> ClusterSequence::exclusive_jets (
const double dcut)
const {
988 int njets = n_exclusive_jets(dcut);
989 return exclusive_jets(njets);
996 vector<PseudoJet> ClusterSequence::exclusive_jets (
const int njets)
const {
1000 if (njets > _initial_n) {
1002 err <<
"Requested " << njets <<
" exclusive jets, but there were only "
1003 << _initial_n <<
" particles in the event";
1004 throw Error(err.str());
1007 return exclusive_jets_up_to(njets);
1013 vector<PseudoJet> ClusterSequence::exclusive_jets_up_to (
const int njets)
const {
1025 (_jet_def.extra_param() <0)) &&
1027 (!_jet_def.plugin()->exclusive_sequence_meaningful()))) {
1028 _exclusive_warnings.warn(
"dcut and exclusive jets for jet-finders other than kt, C/A or genkt with p>=0 should be interpreted with care.");
1035 int stop_point = 2*_initial_n - njets;
1037 if (stop_point < _initial_n) stop_point = _initial_n;
1041 if (2*_initial_n !=
static_cast<int>(_history.size())) {
1043 err <<
"2*_initial_n != _history.size() -- this endangers internal assumptions!\n";
1044 throw Error(err.str());
1053 vector<PseudoJet> jets_local;
1054 for (
unsigned int i = stop_point; i < _history.size(); i++) {
1055 int parent1 = _history[i].parent1;
1056 if (parent1 < stop_point) {
1057 jets_local.push_back(_jets[_history[parent1].jetp_index]);
1059 int parent2 = _history[i].parent2;
1060 if (parent2 < stop_point && parent2 > 0) {
1061 jets_local.push_back(_jets[_history[parent2].jetp_index]);
1067 if (
int(jets_local.size()) != min(_initial_n, njets)) {
1069 err <<
"ClusterSequence::exclusive_jets: size of returned vector ("
1070 <<jets_local.size()<<
") does not coincide with requested number of jets ("
1072 throw Error(err.str());
1081 double ClusterSequence::exclusive_dmerge (
const int njets)
const {
1083 if (njets >= _initial_n) {
return 0.0;}
1084 return _history[2*_initial_n-njets-1].dij;
1093 double ClusterSequence::exclusive_dmerge_max (
const int njets)
const {
1095 if (njets >= _initial_n) {
return 0.0;}
1096 return _history[2*_initial_n-njets-1].max_dij_so_far;
1104 std::vector<PseudoJet> ClusterSequence::exclusive_subjets
1105 (
const PseudoJet & jet,
const double dcut)
const {
1107 set<const history_element*> subhist;
1111 get_subhist_set(subhist, jet, dcut, 0);
1114 vector<PseudoJet> subjets;
1115 subjets.reserve(subhist.size());
1116 for (set<const history_element*>::iterator elem = subhist.begin();
1117 elem != subhist.end(); elem++) {
1118 subjets.push_back(_jets[(*elem)->jetp_index]);
1127 int ClusterSequence::n_exclusive_subjets(
const PseudoJet & jet,
1128 const double dcut)
const {
1129 set<const history_element*> subhist;
1132 get_subhist_set(subhist, jet, dcut, 0);
1133 return subhist.size();
1140 std::vector<PseudoJet> ClusterSequence::exclusive_subjets
1141 (
const PseudoJet & jet,
int nsub)
const {
1142 vector<PseudoJet> subjets = exclusive_subjets_up_to(jet, nsub);
1143 if (
int(subjets.size()) < nsub) {
1145 err <<
"Requested " << nsub <<
" exclusive subjets, but there were only "
1146 << subjets.size() <<
" particles in the jet";
1147 throw Error(err.str());
1157 std::vector<PseudoJet> ClusterSequence::exclusive_subjets_up_to
1158 (
const PseudoJet & jet,
int nsub)
const {
1160 set<const history_element*> subhist;
1163 vector<PseudoJet> subjets;
1164 if (nsub < 0)
throw Error(
"Requested a negative number of subjets. This is nonsensical.");
1165 if (nsub == 0)
return subjets;
1169 get_subhist_set(subhist, jet, -1.0, nsub);
1172 subjets.reserve(subhist.size());
1173 for (set<const history_element*>::iterator elem = subhist.begin();
1174 elem != subhist.end(); elem++) {
1175 subjets.push_back(_jets[(*elem)->jetp_index]);
1186 double ClusterSequence::exclusive_subdmerge(
const PseudoJet & jet,
int nsub)
const {
1187 set<const history_element*> subhist;
1191 get_subhist_set(subhist, jet, -1.0, nsub);
1193 set<const history_element*>::iterator highest = subhist.end();
1197 return (*highest)->dij;
1207 double ClusterSequence::exclusive_subdmerge_max(
const PseudoJet & jet,
int nsub)
const {
1209 set<const history_element*> subhist;
1213 get_subhist_set(subhist, jet, -1.0, nsub);
1215 set<const history_element*>::iterator highest = subhist.end();
1219 return (*highest)->max_dij_so_far;
1231 void ClusterSequence::get_subhist_set(set<const history_element*> & subhist,
1233 double dcut,
int maxjet)
const {
1234 assert(contains(jet));
1244 set<const history_element*>::iterator highest = subhist.end();
1245 assert (highest != subhist.begin());
1249 if (njet == maxjet)
break;
1251 if (elem->parent1 < 0)
break;
1256 subhist.erase(highest);
1257 subhist.insert(&(_history[elem->parent1]));
1258 subhist.insert(&(_history[elem->
parent2]));
1265 bool ClusterSequence::object_in_jet(
const PseudoJet &
object,
1270 assert(contains(
object) && contains(jet));
1272 const PseudoJet * this_object = &object;
1277 }
else if (has_child(*this_object, childp)) {
1278 this_object = childp;
1298 assert ((hist.parent1 >= 0 && hist.
parent2 >= 0) ||
1299 (hist.parent1 < 0 && hist.
parent2 < 0));
1301 if (hist.parent1 < 0) {
1306 parent1 = _jets[_history[hist.parent1].jetp_index];
1307 parent2 = _jets[_history[hist.
parent2].jetp_index];
1309 if (parent1.
perp2() < parent2.
perp2()) std::swap(parent1,parent2);
1329 bool res = has_child(jet, childp);
1346 if (hist.
child >= 0 && _history[hist.
child].jetp_index >= 0) {
1347 childp = &(_jets[_history[hist.
child].jetp_index]);
1367 if (hist.
child >= 0 && _history[hist.
child].parent2 >= 0) {
1372 partner = _jets[_history[child_hist.
parent2].jetp_index];
1375 partner = _jets[_history[child_hist.parent1].jetp_index];
1387 vector<PseudoJet> ClusterSequence::constituents (
const PseudoJet & jet)
const {
1388 vector<PseudoJet> subjets;
1389 add_constituents(jet, subjets);
1402 void ClusterSequence::print_jets_for_root(
const std::vector<PseudoJet> & jets_in,
1403 ostream & ostr)
const {
1404 for (
unsigned i = 0; i < jets_in.size(); i++) {
1406 << jets_in[i].px() <<
" "
1407 << jets_in[i].py() <<
" "
1408 << jets_in[i].pz() <<
" "
1409 << jets_in[i].E() << endl;
1410 vector<PseudoJet> cst = constituents(jets_in[i]);
1411 for (
unsigned j = 0; j < cst.size() ; j++) {
1412 ostr <<
" " << j <<
" "
1413 << cst[j].rap() <<
" "
1414 << cst[j].phi() <<
" "
1415 << cst[j].perp() << endl;
1417 ostr <<
"#END" << endl;
1421 void ClusterSequence::print_jets_for_root(
const std::vector<PseudoJet> & jets_in,
1422 const std::string & filename,
1423 const std::string & comment )
const {
1424 std::ofstream ostr(filename.c_str());
1425 if (comment !=
"") ostr <<
"# " << comment << endl;
1426 print_jets_for_root(jets_in, ostr);
1447 vector<int> ClusterSequence::particle_jet_indices(
1448 const vector<PseudoJet> & jets_in)
const {
1450 vector<int> indices(n_particles());
1453 for (
unsigned ipart = 0; ipart < n_particles(); ipart++)
1454 indices[ipart] = -1;
1458 for (
unsigned ijet = 0; ijet < jets_in.size(); ijet++) {
1460 vector<PseudoJet> jet_constituents(constituents(jets_in[ijet]));
1462 for (
unsigned ip = 0; ip < jet_constituents.size(); ip++) {
1466 unsigned iclust = jet_constituents[ip].cluster_hist_index();
1467 unsigned ipart = history()[iclust].jetp_index;
1468 indices[ipart] = ijet;
1478 void ClusterSequence::add_constituents (
1479 const PseudoJet & jet, vector<PseudoJet> & subjet_vector)
const {
1482 int parent1 = _history[i].parent1;
1483 int parent2 = _history[i].parent2;
1485 if (parent1 == InexistentParent) {
1491 subjet_vector.push_back(_jets[i]);
1496 add_constituents(_jets[_history[parent1].jetp_index], subjet_vector);
1499 if (parent2 != BeamJet) {
1500 add_constituents(_jets[_history[parent2].jetp_index], subjet_vector);
1508 void ClusterSequence::_add_step_to_history (
1511 const int parent2,
const int jetp_index,
1514 history_element element;
1515 element.parent1 = parent1;
1516 element.parent2 = parent2;
1517 element.jetp_index = jetp_index;
1518 element.child = Invalid;
1520 element.max_dij_so_far = max(dij,_history[_history.size()-1].max_dij_so_far);
1521 _history.push_back(element);
1523 int local_step = _history.size()-1;
1535 assert(parent1 >= 0);
1536 if (_history[parent1].child != Invalid){
1537 throw InternalError(
"trying to recomine an object that has previsously been recombined");
1539 _history[parent1].child = local_step;
1541 if (_history[parent2].child != Invalid){
1542 throw InternalError(
"trying to recomine an object that has previsously been recombined");
1544 _history[parent2].child = local_step;
1548 if (jetp_index != Invalid) {
1549 assert(jetp_index >= 0);
1551 _jets[jetp_index].set_cluster_hist_index(local_step);
1552 _set_structure_shared_ptr(_jets[jetp_index]);
1555 if (_writeout_combinations) {
1556 cout << local_step <<
": "
1557 << parent1 <<
" with " << parent2
1558 <<
"; y = "<< dij<<endl;
1571 vector<int> ClusterSequence::unique_history_order()
const {
1577 valarray<int> lowest_constituent(_history.size());
1578 int hist_n = _history.size();
1579 lowest_constituent = hist_n;
1580 for (
int i = 0; i < hist_n; i++) {
1582 lowest_constituent[i] = min(lowest_constituent[i],i);
1584 if (_history[i].child > 0) lowest_constituent[_history[i].child]
1585 = min(lowest_constituent[_history[i].child],lowest_constituent[i]);
1589 valarray<bool> extracted(_history.size()); extracted =
false;
1590 vector<int> unique_tree;
1591 unique_tree.reserve(_history.size());
1594 for (
unsigned i = 0; i < n_particles(); i++) {
1595 if (!extracted[i]) {
1596 unique_tree.push_back(i);
1597 extracted[i] =
true;
1598 _extract_tree_children(i, extracted, lowest_constituent, unique_tree);
1607 void ClusterSequence::_extract_tree_children(
1609 valarray<bool> & extracted,
1610 const valarray<int> & lowest_constituent,
1611 vector<int> & unique_tree)
const {
1612 if (!extracted[position]) {
1615 _extract_tree_parents(position,extracted,lowest_constituent,unique_tree);
1619 int child = _history[position].child;
1620 if (child >= 0) _extract_tree_children(child,extracted,lowest_constituent,unique_tree);
1626 vector<PseudoJet> ClusterSequence::unclustered_particles()
const {
1627 vector<PseudoJet> unclustered;
1628 for (
unsigned i = 0; i < n_particles() ; i++) {
1629 if (_history[i].child == Invalid)
1630 unclustered.push_back(_jets[_history[i].jetp_index]);
1640 vector<PseudoJet> ClusterSequence::childless_pseudojets()
const {
1641 vector<PseudoJet> unclustered;
1642 for (
unsigned i = 0; i < _history.size() ; i++) {
1643 if ((_history[i].child == Invalid) && (_history[i].parent2 != BeamJet))
1644 unclustered.push_back(_jets[_history[i].jetp_index]);
1655 bool ClusterSequence::contains(
const PseudoJet & jet)
const {
1666 void ClusterSequence::_extract_tree_parents(
1668 valarray<bool> & extracted,
1669 const valarray<int> & lowest_constituent,
1670 vector<int> & unique_tree)
const {
1672 if (!extracted[position]) {
1673 int parent1 = _history[position].parent1;
1674 int parent2 = _history[position].parent2;
1677 if (parent1 >= 0 && parent2 >= 0) {
1678 if (lowest_constituent[parent1] > lowest_constituent[parent2])
1679 std::swap(parent1, parent2);
1682 if (parent1 >= 0 && !extracted[parent1])
1683 _extract_tree_parents(parent1,extracted,lowest_constituent,unique_tree);
1684 if (parent2 >= 0 && !extracted[parent2])
1685 _extract_tree_parents(parent2,extracted,lowest_constituent,unique_tree);
1688 unique_tree.push_back(position);
1689 extracted[position] =
true;
1698 void ClusterSequence::_do_ij_recombination_step(
1699 const int jet_i,
const int jet_j,
1709 _jet_def.recombiner()->recombine(_jets[jet_i], _jets[jet_j], newjet);
1710 _jets.push_back(newjet);
1715 newjet_k = _jets.size()-1;
1718 int newstep_k = _history.size();
1720 _jets[newjet_k].set_cluster_hist_index(newstep_k);
1723 int hist_i = _jets[jet_i].cluster_hist_index();
1724 int hist_j = _jets[jet_j].cluster_hist_index();
1726 _add_step_to_history(min(hist_i, hist_j), max(hist_i,hist_j),
1739 void ClusterSequence::_do_iB_recombination_step(
1740 const int jet_i,
const double diB) {
1742 _add_step_to_history(_jets[jet_i].cluster_hist_index(),BeamJet,
1759 void ClusterSequence::_set_structure_shared_ptr(
PseudoJet & j) {
1763 _update_structure_use_count();
1768 void ClusterSequence::_update_structure_use_count() {
1771 _structure_use_count_after_construction = _structure_shared_ptr.use_count();
1778 void ClusterSequence::delete_self_when_unused() {
1786 int new_count = _structure_shared_ptr.use_count() - _structure_use_count_after_construction;
1787 if (new_count <= 0) {
1788 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");
1795 _structure_shared_ptr.set_count(new_count);
1797 _deletes_self_when_unused =
true;
1801 FASTJET_END_NAMESPACE