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"
39 #include "fastjet/internal/LazyTiling9SeparateGhosts.hh"
49 FASTJET_BEGIN_NAMESPACE
145 std::ostream * ClusterSequence::_fastjet_banner_ostr = &cout;
149 ClusterSequence::~ClusterSequence () {
152 if (_structure_shared_ptr()){
153 ClusterSequenceStructure* csi =
dynamic_cast<ClusterSequenceStructure*
>(_structure_shared_ptr());
159 csi->set_associated_cs(NULL);
167 if (_deletes_self_when_unused) {
168 _structure_shared_ptr.set_count(_structure_shared_ptr.use_count()
169 + _structure_use_count_after_construction);
175 void ClusterSequence::signal_imminent_self_deletion()
const {
189 assert(_deletes_self_when_unused);
190 _deletes_self_when_unused =
false;
205 void ClusterSequence::_initialise_and_run (
207 const bool & writeout_combinations) {
210 _decant_options(jet_def_in, writeout_combinations);
213 _initialise_and_run_no_decant();
217 void ClusterSequence::_initialise_and_run_no_decant () {
221 _fill_initial_history();
224 if (n_particles() == 0)
return;
227 if (_jet_algorithm == plugin_algorithm) {
229 _plugin_activated =
true;
231 _jet_def.plugin()->run_clustering( (*
this) );
232 _plugin_activated =
false;
233 _update_structure_use_count();
235 }
else if (_jet_algorithm == ee_kt_algorithm ||
236 _jet_algorithm == ee_genkt_algorithm) {
239 if (_jet_algorithm == ee_kt_algorithm) {
243 assert(_Rparam > 2.0);
257 _R2 = 2 * ( 3.0 + cos(_Rparam) );
259 _R2 = 2 * ( 1.0 - cos(_Rparam) );
263 _simple_N2_cluster_EEBriefJet();
265 }
else if (_jet_algorithm == undefined_jet_algorithm) {
266 throw Error(
"A ClusterSequence cannot be created with an uninitialised JetDefinition");
280 if (_strategy == Best) {
281 _strategy = _best_strategy();
286 }
else if (_strategy == BestFJ30) {
287 int N = _jets.size();
293 if (min(1.0,max(0.1,_Rparam)*3.3)*N <= 30) {
298 }
else if ((N > 16000/pow(_Rparam,1.15) && _jet_def.jet_algorithm() !=
antikt_algorithm)
299 || N > 35000/pow(_Rparam,1.15)) {
302 }
else if (N <= 450) {
313 if (_Rparam >= twopi) {
314 if ( _strategy == NlnN
315 || _strategy == NlnN3pi
316 || _strategy == NlnNCam
317 || _strategy == NlnNCam2pi2R
318 || _strategy == NlnNCam4pi) {
325 if (_jet_def.strategy() !=
Best && _strategy != _jet_def.strategy()) {
327 oss <<
"Cluster strategy " << strategy_string(_jet_def.strategy())
328 <<
" automatically changed to " << strategy_string()
329 <<
" because the former is not supported for R = " << _Rparam
331 _changed_strategy_warning.warn(oss.str());
341 if (_strategy == N2Plain) {
343 this->_simple_N2_cluster_BriefJet();
344 }
else if (_strategy == N2Tiled) {
345 this->_faster_tiled_N2_cluster();
346 }
else if (_strategy == N2MinHeapTiled) {
347 this->_minheap_faster_tiled_N2_cluster();
348 }
else if (_strategy == N2MHTLazy9Alt) {
351 _plugin_activated =
true;
352 LazyTiling9Alt tiling(*
this);
354 _plugin_activated =
false;
356 }
else if (_strategy == N2MHTLazy25) {
359 _plugin_activated =
true;
360 LazyTiling25 tiling(*
this);
362 _plugin_activated =
false;
364 }
else if (_strategy == N2MHTLazy9) {
367 _plugin_activated =
true;
368 LazyTiling9 tiling(*
this);
370 _plugin_activated =
false;
372 }
else if (_strategy == N2MHTLazy9AntiKtSeparateGhosts) {
375 _plugin_activated =
true;
376 LazyTiling9SeparateGhosts tiling(*
this);
378 _plugin_activated =
false;
380 }
else if (_strategy == NlnN) {
381 this->_delaunay_cluster();
382 }
else if (_strategy == NlnNCam) {
383 this->_CP2DChan_cluster_2piMultD();
384 }
else if (_strategy == NlnN3pi || _strategy == NlnN4pi ) {
385 this->_delaunay_cluster();
386 }
else if (_strategy == N3Dumb ) {
387 this->_really_dumb_cluster();
388 }
else if (_strategy == N2PoorTiled) {
389 this->_tiled_N2_cluster();
390 }
else if (_strategy == NlnNCam4pi) {
391 this->_CP2DChan_cluster();
392 }
else if (_strategy == NlnNCam2pi2R) {
393 this->_CP2DChan_cluster_2pi2R();
396 err <<
"Unrecognised value for strategy: "<<_strategy;
397 throw Error(err.str());
404 bool ClusterSequence::_first_time =
true;
405 LimitedWarning ClusterSequence::_exclusive_warnings;
411 return "FastJet version "+string(fastjet_version);
417 void ClusterSequence::print_banner() {
419 if (!_first_time) {
return;}
423 ostream * ostr = _fastjet_banner_ostr;
426 (*ostr) <<
"#--------------------------------------------------------------------------\n";
427 (*ostr) <<
"# FastJet release " << fastjet_version << endl;
428 (*ostr) <<
"# M. Cacciari, G.P. Salam and G. Soyez \n";
429 (*ostr) <<
"# A software package for jet finding and analysis at colliders \n";
430 (*ostr) <<
"# http://fastjet.fr \n";
432 (*ostr) <<
"# Please cite EPJC72(2012)1896 [arXiv:1111.6097] if you use this package\n";
433 (*ostr) <<
"# for scientific work and optionally PLB641(2006)57 [hep-ph/0512210]. \n";
435 (*ostr) <<
"# FastJet is provided without warranty under the terms of the GNU GPLv2.\n";
436 (*ostr) <<
"# It uses T. Chan's closest pair algorithm, S. Fortune's Voronoi code";
438 (*ostr) <<
",\n# CGAL ";
442 (*ostr) <<
"and 3rd party plugin jet algorithms. See COPYING file for details.\n";
443 (*ostr) <<
"#--------------------------------------------------------------------------\n";
451 const bool & writeout_combinations) {
453 _jet_def = jet_def_in;
454 _writeout_combinations = writeout_combinations;
458 _decant_options_partial();
463 void ClusterSequence::_decant_options_partial() {
467 _jet_algorithm = _jet_def.jet_algorithm();
468 _Rparam = _jet_def.R(); _R2 = _Rparam*_Rparam; _invR2 = 1.0/_R2;
469 _strategy = _jet_def.strategy();
472 _plugin_activated =
false;
476 _update_structure_use_count();
482 void ClusterSequence::_fill_initial_history () {
487 _jets.reserve(_jets.size()*2);
488 _history.reserve(_jets.size()*2);
492 for (
int i = 0; i < static_cast<int>(_jets.size()) ; i++) {
494 element.parent1 = InexistentParent;
495 element.
parent2 = InexistentParent;
496 element.
child = Invalid;
501 _history.push_back(element);
504 _jet_def.recombiner()->preprocess(_jets[i]);
507 _jets[i].set_cluster_hist_index(i);
508 _set_structure_shared_ptr(_jets[i]);
511 _Qtot += _jets[i].E();
513 _initial_n = _jets.size();
514 _deletes_self_when_unused =
false;
519 string ClusterSequence::strategy_string (
Strategy strategy_in)
const {
521 switch(strategy_in) {
523 strategy =
"NlnN";
break;
525 strategy =
"NlnN3pi";
break;
527 strategy =
"NlnN4pi";
break;
529 strategy =
"N2Plain";
break;
531 strategy =
"N2Tiled";
break;
533 strategy =
"N2MinHeapTiled";
break;
535 strategy =
"N2PoorTiled";
break;
537 strategy =
"N2MHTLazy9";
break;
539 strategy =
"N2MHTLazy9Alt";
break;
541 strategy =
"N2MHTLazy25";
break;
543 strategy =
"N2MHTLazy9AntiKtSeparateGhosts";
break;
545 strategy =
"N3Dumb";
break;
547 strategy =
"NlnNCam4pi";
break;
549 strategy =
"NlnNCam2pi2R";
break;
551 strategy =
"NlnNCam";
break;
553 strategy =
"plugin strategy";
break;
555 strategy =
"Unrecognized";
561 double ClusterSequence::jet_scale_for_algorithm(
566 double kt2=jet.
kt2();
567 return kt2 > 1e-300 ? 1.0/kt2 : 1e300;
569 double kt2 = jet.
kt2();
570 double p = jet_def().extra_param();
571 if (p <= 0 && kt2 < 1e-300) kt2 = 1e-300;
574 double kt2 = jet.
kt2();
575 double lim = _jet_def.extra_param();
576 if (kt2 < lim*lim && kt2 != 0.0) {
579 }
else {
throw Error(
"Unrecognised jet algorithm");}
598 int N = _jets.size();
601 double bounded_R = max(_Rparam, 0.1);
605 if (N <= 30 || N <= 39.0/(bounded_R + 0.6)) {
616 const static _Parabola N_Tiled_to_MHT_lowR (-45.4947,54.3528,44.6283);
617 const static _Parabola L_MHT_to_MHTLazy9_lowR (0.677807,-1.05006,10.6994);
618 const static _Parabola L_MHTLazy9_to_MHTLazy25_akt_lowR(0.169967,-0.512589,12.1572);
619 const static _Parabola L_MHTLazy9_to_MHTLazy25_kt_lowR (0.16237,-0.484612,12.3373);
620 const static _Parabola L_MHTLazy9_to_MHTLazy25_cam_lowR = L_MHTLazy9_to_MHTLazy25_kt_lowR;
621 const static _Parabola L_MHTLazy25_to_NlnN_akt_lowR (0.0472051,-0.22043,15.9196);
622 const static _Parabola L_MHTLazy25_to_NlnN_kt_lowR (0.118609,-0.326811,14.8287);
623 const static _Parabola L_MHTLazy25_to_NlnN_cam_lowR (0.10119,-0.295748,14.3924);
625 const static _Line L_Tiled_to_MHTLazy9_medR (-1.31304,7.29621);
626 const static _Parabola L_MHTLazy9_to_MHTLazy25_akt_medR = L_MHTLazy9_to_MHTLazy25_akt_lowR;
627 const static _Parabola L_MHTLazy9_to_MHTLazy25_kt_medR = L_MHTLazy9_to_MHTLazy25_kt_lowR;
628 const static _Parabola L_MHTLazy9_to_MHTLazy25_cam_medR = L_MHTLazy9_to_MHTLazy25_cam_lowR;
629 const static _Parabola L_MHTLazy25_to_NlnN_akt_medR = L_MHTLazy25_to_NlnN_akt_lowR;
630 const static _Parabola L_MHTLazy25_to_NlnN_kt_medR = L_MHTLazy25_to_NlnN_kt_lowR;
631 const static _Parabola L_MHTLazy25_to_NlnN_cam_medR = L_MHTLazy25_to_NlnN_cam_lowR;
633 const static double N_Plain_to_MHTLazy9_largeR = 75;
634 const static double N_MHTLazy9_to_MHTLazy25_akt_largeR = 700;
635 const static double N_MHTLazy9_to_MHTLazy25_kt_largeR = 1000;
636 const static double N_MHTLazy9_to_MHTLazy25_cam_largeR = 1000;
637 const static double N_MHTLazy25_to_NlnN_akt_largeR = 100000;
638 const static double N_MHTLazy25_to_NlnN_kt_largeR = 40000;
639 const static double N_MHTLazy25_to_NlnN_cam_largeR = 15000;
650 double p = jet_def().extra_param();
658 jet_algorithm = _jet_algorithm;
661 if (bounded_R < 0.65) {
663 if (N < N_Tiled_to_MHT_lowR(bounded_R))
return N2Tiled;
664 double logN = log(
double(N));
665 if (logN < L_MHT_to_MHTLazy9_lowR(bounded_R))
return N2MinHeapTiled;
668 if (logN < L_MHTLazy9_to_MHTLazy25_akt_lowR(bounded_R))
return N2MHTLazy9;
669 else if (logN < L_MHTLazy25_to_NlnN_akt_lowR(bounded_R))
return N2MHTLazy25;
672 if (logN < L_MHTLazy9_to_MHTLazy25_kt_lowR(bounded_R))
return N2MHTLazy9;
673 else if (logN < L_MHTLazy25_to_NlnN_kt_lowR(bounded_R))
return N2MHTLazy25;
676 if (logN < L_MHTLazy9_to_MHTLazy25_cam_lowR(bounded_R))
return N2MHTLazy9;
677 else if (logN < L_MHTLazy25_to_NlnN_cam_lowR(bounded_R))
return N2MHTLazy25;
681 }
else if (bounded_R < 0.5*pi) {
683 double logN = log(
double(N));
684 if (logN < L_Tiled_to_MHTLazy9_medR(bounded_R))
return N2Tiled;
687 if (logN < L_MHTLazy9_to_MHTLazy25_akt_medR(bounded_R))
return N2MHTLazy9;
688 else if (logN < L_MHTLazy25_to_NlnN_akt_medR(bounded_R))
return N2MHTLazy25;
691 if (logN < L_MHTLazy9_to_MHTLazy25_kt_medR(bounded_R))
return N2MHTLazy9;
692 else if (logN < L_MHTLazy25_to_NlnN_kt_medR(bounded_R))
return N2MHTLazy25;
695 if (logN < L_MHTLazy9_to_MHTLazy25_cam_medR(bounded_R))
return N2MHTLazy9;
696 else if (logN < L_MHTLazy25_to_NlnN_cam_medR(bounded_R))
return N2MHTLazy25;
702 if (N < N_Plain_to_MHTLazy9_largeR)
return N2Plain;
705 if (N < N_MHTLazy9_to_MHTLazy25_akt_largeR)
return N2MHTLazy9;
706 else if (N < N_MHTLazy25_to_NlnN_akt_largeR)
return N2MHTLazy25;
709 if (N < N_MHTLazy9_to_MHTLazy25_kt_largeR)
return N2MHTLazy9;
710 else if (N < N_MHTLazy25_to_NlnN_kt_largeR)
return N2MHTLazy25;
713 if (N < N_MHTLazy9_to_MHTLazy25_cam_largeR)
return N2MHTLazy9;
714 else if (N < N_MHTLazy25_to_NlnN_cam_largeR)
return N2MHTLazy25;
720 bool code_should_never_reach_here =
false;
721 assert(code_should_never_reach_here);
784 if (will_delete_self_when_unused())
785 throw(
Error(
"cannot use CS::transfer_from_sequence after a call to delete_self_when_unused()"));
788 _jet_def = from_seq._jet_def ;
789 _writeout_combinations = from_seq._writeout_combinations ;
790 _initial_n = from_seq._initial_n ;
791 _Rparam = from_seq._Rparam ;
793 _invR2 = from_seq._invR2 ;
794 _strategy = from_seq._strategy ;
795 _jet_algorithm = from_seq._jet_algorithm ;
796 _plugin_activated = from_seq._plugin_activated ;
802 _jets = (*action_on_jets)(from_seq.
_jets);
804 _jets = from_seq.
_jets;
808 _extras = from_seq._extras;
811 if (_structure_shared_ptr()) {
815 if (_deletes_self_when_unused)
throw Error(
"transfer_from_sequence cannot be used for a cluster sequence that deletes self when unused");
826 _update_structure_use_count();
828 for (
unsigned int i=0; i<_jets.size(); i++){
831 _jets[i].set_cluster_hist_index(from_seq.
_jets[i].cluster_hist_index());
834 _set_structure_shared_ptr(_jets[i]);
842 void ClusterSequence::plugin_record_ij_recombination(
843 int jet_i,
int jet_j,
double dij,
844 const PseudoJet & newjet,
int & newjet_k) {
846 plugin_record_ij_recombination(jet_i, jet_j, dij, newjet_k);
849 int tmp_index = _jets[newjet_k].cluster_hist_index();
850 _jets[newjet_k] = newjet;
852 _set_structure_shared_ptr(_jets[newjet_k]);
858 vector<PseudoJet> ClusterSequence::inclusive_jets (
const double ptmin)
const{
859 double dcut = ptmin*ptmin;
860 int i = _history.size() - 1;
861 vector<PseudoJet> jets_local;
867 if (_history[i].max_dij_so_far < dcut) {
break;}
868 if (_history[i].parent2 == BeamJet && _history[i].dij >= dcut) {
870 int parent1 = _history[i].parent1;
871 jets_local.push_back(_jets[_history[parent1].jetp_index]);}
879 if (_history[i].parent2 != BeamJet) {
break;}
880 int parent1 = _history[i].parent1;
881 const PseudoJet & jet = _jets[_history[parent1].jetp_index];
882 if (jet.
perp2() >= dcut) {jets_local.push_back(jet);}
895 if (_history[i].parent2 == BeamJet) {
896 int parent1 = _history[i].parent1;
897 const PseudoJet & jet = _jets[_history[parent1].jetp_index];
898 if (jet.
perp2() >= dcut) {jets_local.push_back(jet);}
902 }
else {
throw Error(
"cs::inclusive_jets(...): Unrecognized jet algorithm");}
910 int ClusterSequence::n_exclusive_jets (
const double dcut)
const {
914 int i = _history.size() - 1;
916 if (_history[i].max_dij_so_far <= dcut) {
break;}
919 int stop_point = i + 1;
922 int njets = 2*_initial_n - stop_point;
929 vector<PseudoJet> ClusterSequence::exclusive_jets (
const double dcut)
const {
930 int njets = n_exclusive_jets(dcut);
931 return exclusive_jets(njets);
938 vector<PseudoJet> ClusterSequence::exclusive_jets (
const int njets)
const {
942 if (njets > _initial_n) {
944 err <<
"Requested " << njets <<
" exclusive jets, but there were only "
945 << _initial_n <<
" particles in the event";
946 throw Error(err.str());
949 return exclusive_jets_up_to(njets);
955 vector<PseudoJet> ClusterSequence::exclusive_jets_up_to (
const int njets)
const {
967 (_jet_def.extra_param() <0)) &&
969 (!_jet_def.plugin()->exclusive_sequence_meaningful()))) {
970 _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.");
977 int stop_point = 2*_initial_n - njets;
979 if (stop_point < _initial_n) stop_point = _initial_n;
983 if (2*_initial_n != static_cast<int>(_history.size())) {
985 err <<
"2*_initial_n != _history.size() -- this endangers internal assumptions!\n";
986 throw Error(err.str());
995 vector<PseudoJet> jets_local;
996 for (
unsigned int i = stop_point; i < _history.size(); i++) {
997 int parent1 = _history[i].parent1;
998 if (parent1 < stop_point) {
999 jets_local.push_back(_jets[_history[parent1].jetp_index]);
1001 int parent2 = _history[i].parent2;
1002 if (parent2 < stop_point && parent2 > 0) {
1003 jets_local.push_back(_jets[_history[parent2].jetp_index]);
1009 if (
int(jets_local.size()) != min(_initial_n, njets)) {
1011 err <<
"ClusterSequence::exclusive_jets: size of returned vector ("
1012 <<jets_local.size()<<
") does not coincide with requested number of jets ("
1014 throw Error(err.str());
1023 double ClusterSequence::exclusive_dmerge (
const int njets)
const {
1025 if (njets >= _initial_n) {
return 0.0;}
1026 return _history[2*_initial_n-njets-1].dij;
1035 double ClusterSequence::exclusive_dmerge_max (
const int njets)
const {
1037 if (njets >= _initial_n) {
return 0.0;}
1038 return _history[2*_initial_n-njets-1].max_dij_so_far;
1046 std::vector<PseudoJet> ClusterSequence::exclusive_subjets
1049 set<const history_element*> subhist;
1053 get_subhist_set(subhist, jet, dcut, 0);
1056 vector<PseudoJet> subjets;
1057 subjets.reserve(subhist.size());
1058 for (set<const history_element*>::iterator elem = subhist.begin();
1059 elem != subhist.end(); elem++) {
1060 subjets.push_back(_jets[(*elem)->jetp_index]);
1069 int ClusterSequence::n_exclusive_subjets(
const PseudoJet & jet,
1070 const double dcut)
const {
1071 set<const history_element*> subhist;
1074 get_subhist_set(subhist, jet, dcut, 0);
1075 return subhist.size();
1082 std::vector<PseudoJet> ClusterSequence::exclusive_subjets
1084 vector<PseudoJet> subjets = exclusive_subjets_up_to(jet, nsub);
1085 if (
int(subjets.size()) < nsub) {
1087 err <<
"Requested " << nsub <<
" exclusive subjets, but there were only "
1088 << subjets.size() <<
" particles in the jet";
1089 throw Error(err.str());
1099 std::vector<PseudoJet> ClusterSequence::exclusive_subjets_up_to
1102 set<const history_element*> subhist;
1105 vector<PseudoJet> subjets;
1106 if (nsub < 0)
throw Error(
"Requested a negative number of subjets. This is nonsensical.");
1107 if (nsub == 0)
return subjets;
1111 get_subhist_set(subhist, jet, -1.0, nsub);
1114 subjets.reserve(subhist.size());
1115 for (set<const history_element*>::iterator elem = subhist.begin();
1116 elem != subhist.end(); elem++) {
1117 subjets.push_back(_jets[(*elem)->jetp_index]);
1128 double ClusterSequence::exclusive_subdmerge(
const PseudoJet & jet,
int nsub)
const {
1129 set<const history_element*> subhist;
1133 get_subhist_set(subhist, jet, -1.0, nsub);
1135 set<const history_element*>::iterator highest = subhist.end();
1139 return (*highest)->dij;
1149 double ClusterSequence::exclusive_subdmerge_max(
const PseudoJet & jet,
int nsub)
const {
1151 set<const history_element*> subhist;
1155 get_subhist_set(subhist, jet, -1.0, nsub);
1157 set<const history_element*>::iterator highest = subhist.end();
1161 return (*highest)->max_dij_so_far;
1173 void ClusterSequence::get_subhist_set(set<const history_element*> & subhist,
1175 double dcut,
int maxjet)
const {
1176 assert(contains(jet));
1186 set<const history_element*>::iterator highest = subhist.end();
1187 assert (highest != subhist.begin());
1191 if (njet == maxjet)
break;
1193 if (elem->parent1 < 0)
break;
1198 subhist.erase(highest);
1199 subhist.insert(&(_history[elem->parent1]));
1200 subhist.insert(&(_history[elem->
parent2]));
1207 bool ClusterSequence::object_in_jet(
const PseudoJet &
object,
1212 assert(contains(
object) && contains(jet));
1214 const PseudoJet * this_object = &object;
1219 }
else if (has_child(*this_object, childp)) {
1220 this_object = childp;
1240 assert ((hist.parent1 >= 0 && hist.
parent2 >= 0) ||
1241 (hist.parent1 < 0 && hist.
parent2 < 0));
1243 if (hist.parent1 < 0) {
1248 parent1 = _jets[_history[hist.parent1].jetp_index];
1249 parent2 = _jets[_history[hist.
parent2].jetp_index];
1251 if (parent1.
perp2() < parent2.
perp2()) std::swap(parent1,parent2);
1271 bool res = has_child(jet, childp);
1288 if (hist.
child >= 0 && _history[hist.
child].jetp_index >= 0) {
1289 childp = &(_jets[_history[hist.
child].jetp_index]);
1309 if (hist.
child >= 0 && _history[hist.
child].parent2 >= 0) {
1314 partner = _jets[_history[child_hist.
parent2].jetp_index];
1317 partner = _jets[_history[child_hist.parent1].jetp_index];
1329 vector<PseudoJet> ClusterSequence::constituents (
const PseudoJet & jet)
const {
1330 vector<PseudoJet> subjets;
1331 add_constituents(jet, subjets);
1344 void ClusterSequence::print_jets_for_root(
const std::vector<PseudoJet> & jets_in,
1345 ostream & ostr)
const {
1346 for (
unsigned i = 0; i < jets_in.size(); i++) {
1348 << jets_in[i].px() <<
" "
1349 << jets_in[i].py() <<
" "
1350 << jets_in[i].pz() <<
" "
1351 << jets_in[i].E() << endl;
1352 vector<PseudoJet> cst = constituents(jets_in[i]);
1353 for (
unsigned j = 0; j < cst.size() ; j++) {
1354 ostr <<
" " << j <<
" "
1355 << cst[j].rap() <<
" "
1356 << cst[j].phi() <<
" "
1357 << cst[j].perp() << endl;
1359 ostr <<
"#END" << endl;
1363 void ClusterSequence::print_jets_for_root(
const std::vector<PseudoJet> & jets_in,
1364 const std::string & filename,
1365 const std::string & comment )
const {
1366 std::ofstream ostr(filename.c_str());
1367 if (comment !=
"") ostr <<
"# " << comment << endl;
1368 print_jets_for_root(jets_in, ostr);
1389 vector<int> ClusterSequence::particle_jet_indices(
1390 const vector<PseudoJet> & jets_in)
const {
1392 vector<int> indices(n_particles());
1395 for (
unsigned ipart = 0; ipart < n_particles(); ipart++)
1396 indices[ipart] = -1;
1400 for (
unsigned ijet = 0; ijet < jets_in.size(); ijet++) {
1402 vector<PseudoJet> jet_constituents(constituents(jets_in[ijet]));
1404 for (
unsigned ip = 0; ip < jet_constituents.size(); ip++) {
1408 unsigned iclust = jet_constituents[ip].cluster_hist_index();
1409 unsigned ipart = history()[iclust].jetp_index;
1410 indices[ipart] = ijet;
1420 void ClusterSequence::add_constituents (
1421 const PseudoJet & jet, vector<PseudoJet> & subjet_vector)
const {
1424 int parent1 = _history[i].parent1;
1425 int parent2 = _history[i].parent2;
1427 if (parent1 == InexistentParent) {
1433 subjet_vector.push_back(_jets[i]);
1438 add_constituents(_jets[_history[parent1].jetp_index], subjet_vector);
1441 if (parent2 != BeamJet) {
1442 add_constituents(_jets[_history[parent2].jetp_index], subjet_vector);
1450 void ClusterSequence::_add_step_to_history (
1451 const int step_number,
const int parent1,
1452 const int parent2,
const int jetp_index,
1455 history_element element;
1456 element.parent1 = parent1;
1457 element.parent2 = parent2;
1458 element.jetp_index = jetp_index;
1459 element.child = Invalid;
1461 element.max_dij_so_far = max(dij,_history[_history.size()-1].max_dij_so_far);
1462 _history.push_back(element);
1464 int local_step = _history.size()-1;
1465 assert(local_step == step_number);
1467 assert(parent1 >= 0);
1468 _history[parent1].child = local_step;
1469 if (parent2 >= 0) {_history[parent2].child = local_step;}
1472 if (jetp_index != Invalid) {
1473 assert(jetp_index >= 0);
1475 _jets[jetp_index].set_cluster_hist_index(local_step);
1476 _set_structure_shared_ptr(_jets[jetp_index]);
1479 if (_writeout_combinations) {
1480 cout << local_step <<
": "
1481 << parent1 <<
" with " << parent2
1482 <<
"; y = "<< dij<<endl;
1495 vector<int> ClusterSequence::unique_history_order()
const {
1501 valarray<int> lowest_constituent(_history.size());
1502 int hist_n = _history.size();
1503 lowest_constituent = hist_n;
1504 for (
int i = 0; i < hist_n; i++) {
1506 lowest_constituent[i] = min(lowest_constituent[i],i);
1508 if (_history[i].child > 0) lowest_constituent[_history[i].child]
1509 = min(lowest_constituent[_history[i].child],lowest_constituent[i]);
1513 valarray<bool> extracted(_history.size()); extracted =
false;
1514 vector<int> unique_tree;
1515 unique_tree.reserve(_history.size());
1518 for (
unsigned i = 0; i < n_particles(); i++) {
1519 if (!extracted[i]) {
1520 unique_tree.push_back(i);
1521 extracted[i] =
true;
1522 _extract_tree_children(i, extracted, lowest_constituent, unique_tree);
1531 void ClusterSequence::_extract_tree_children(
1533 valarray<bool> & extracted,
1534 const valarray<int> & lowest_constituent,
1535 vector<int> & unique_tree)
const {
1536 if (!extracted[position]) {
1539 _extract_tree_parents(position,extracted,lowest_constituent,unique_tree);
1543 int child = _history[position].child;
1544 if (child >= 0) _extract_tree_children(child,extracted,lowest_constituent,unique_tree);
1550 vector<PseudoJet> ClusterSequence::unclustered_particles()
const {
1551 vector<PseudoJet> unclustered;
1552 for (
unsigned i = 0; i < n_particles() ; i++) {
1553 if (_history[i].child == Invalid)
1554 unclustered.push_back(_jets[_history[i].jetp_index]);
1564 vector<PseudoJet> ClusterSequence::childless_pseudojets()
const {
1565 vector<PseudoJet> unclustered;
1566 for (
unsigned i = 0; i < _history.size() ; i++) {
1567 if ((_history[i].child == Invalid) && (_history[i].parent2 != BeamJet))
1568 unclustered.push_back(_jets[_history[i].jetp_index]);
1579 bool ClusterSequence::contains(
const PseudoJet & jet)
const {
1590 void ClusterSequence::_extract_tree_parents(
1592 valarray<bool> & extracted,
1593 const valarray<int> & lowest_constituent,
1594 vector<int> & unique_tree)
const {
1596 if (!extracted[position]) {
1597 int parent1 = _history[position].parent1;
1598 int parent2 = _history[position].parent2;
1601 if (parent1 >= 0 && parent2 >= 0) {
1602 if (lowest_constituent[parent1] > lowest_constituent[parent2])
1603 std::swap(parent1, parent2);
1606 if (parent1 >= 0 && !extracted[parent1])
1607 _extract_tree_parents(parent1,extracted,lowest_constituent,unique_tree);
1608 if (parent2 >= 0 && !extracted[parent2])
1609 _extract_tree_parents(parent2,extracted,lowest_constituent,unique_tree);
1612 unique_tree.push_back(position);
1613 extracted[position] =
true;
1622 void ClusterSequence::_do_ij_recombination_step(
1623 const int jet_i,
const int jet_j,
1633 _jet_def.recombiner()->recombine(_jets[jet_i], _jets[jet_j], newjet);
1634 _jets.push_back(newjet);
1639 newjet_k = _jets.size()-1;
1642 int newstep_k = _history.size();
1644 _jets[newjet_k].set_cluster_hist_index(newstep_k);
1647 int hist_i = _jets[jet_i].cluster_hist_index();
1648 int hist_j = _jets[jet_j].cluster_hist_index();
1650 _add_step_to_history(newstep_k, min(hist_i, hist_j), max(hist_i,hist_j),
1659 void ClusterSequence::_do_iB_recombination_step(
1660 const int jet_i,
const double diB) {
1662 int newstep_k = _history.size();
1665 _add_step_to_history(newstep_k,_jets[jet_i].cluster_hist_index(),BeamJet,
1676 void ClusterSequence::_set_structure_shared_ptr(
PseudoJet & j) {
1680 _update_structure_use_count();
1685 void ClusterSequence::_update_structure_use_count() {
1688 _structure_use_count_after_construction = _structure_shared_ptr.use_count();
1695 void ClusterSequence::delete_self_when_unused() {
1703 int new_count = _structure_shared_ptr.use_count() - _structure_use_count_after_construction;
1704 if (new_count <= 0) {
1705 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");
1708 _structure_shared_ptr.set_count(new_count);
1709 _deletes_self_when_unused =
true;
1713 FASTJET_END_NAMESPACE