29 #include "fastjet/Error.hh"
30 #include "fastjet/PseudoJet.hh"
31 #include "fastjet/ClusterSequence.hh"
32 #include "fastjet/ClusterSequenceStructure.hh"
33 #include "fastjet/version.hh"
43 FASTJET_BEGIN_NAMESPACE
139 std::ostream * ClusterSequence::_fastjet_banner_ostr = &cout;
143 ClusterSequence::~ClusterSequence () {
146 if (_structure_shared_ptr()){
147 ClusterSequenceStructure* csi =
dynamic_cast<ClusterSequenceStructure*
>(_structure_shared_ptr());
153 csi->set_associated_cs(NULL);
161 if (_deletes_self_when_unused) {
162 _structure_shared_ptr.set_count(_structure_shared_ptr.use_count()
163 + _structure_use_count_after_construction);
169 void ClusterSequence::signal_imminent_self_deletion()
const {
183 assert(_deletes_self_when_unused);
184 _deletes_self_when_unused =
false;
199 void ClusterSequence::_initialise_and_run (
201 const bool & writeout_combinations) {
204 _decant_options(jet_def_in, writeout_combinations);
207 _initialise_and_run_no_decant();
211 void ClusterSequence::_initialise_and_run_no_decant () {
215 _fill_initial_history();
218 if (n_particles() == 0)
return;
221 if (_jet_algorithm == plugin_algorithm) {
223 _plugin_activated =
true;
225 _jet_def.plugin()->run_clustering( (*
this) );
226 _plugin_activated =
false;
227 _update_structure_use_count();
229 }
else if (_jet_algorithm == ee_kt_algorithm ||
230 _jet_algorithm == ee_genkt_algorithm) {
233 if (_jet_algorithm == ee_kt_algorithm) {
237 assert(_Rparam > 2.0);
251 _R2 = 2 * ( 3.0 + cos(_Rparam) );
253 _R2 = 2 * ( 1.0 - cos(_Rparam) );
257 _simple_N2_cluster_EEBriefJet();
259 }
else if (_jet_algorithm == undefined_jet_algorithm) {
260 throw Error(
"A ClusterSequence cannot be created with an uninitialised JetDefinition");
274 if (_strategy == Best) {
275 int N = _jets.size();
281 if (min(1.0,max(0.1,_Rparam)*3.3)*N <= 30) {
286 }
else if ((N > 16000/pow(_Rparam,1.15) && _jet_def.jet_algorithm() !=
antikt_algorithm)
287 || N > 35000/pow(_Rparam,1.15)) {
290 }
else if (N <= 450) {
301 if (_Rparam >= twopi) {
302 if ( _strategy == NlnN
303 || _strategy == NlnN3pi
304 || _strategy == NlnNCam
305 || _strategy == NlnNCam2pi2R
306 || _strategy == NlnNCam4pi) {
313 if (_jet_def.strategy() !=
Best && _strategy != _jet_def.strategy()) {
315 oss <<
"Cluster strategy " << strategy_string(_jet_def.strategy())
316 <<
" automatically changed to " << strategy_string()
317 <<
" because the former is not supported for R = " << _Rparam
319 _changed_strategy_warning.warn(oss.str());
329 if (_strategy == N2Plain) {
331 this->_simple_N2_cluster_BriefJet();
332 }
else if (_strategy == N2Tiled) {
333 this->_faster_tiled_N2_cluster();
334 }
else if (_strategy == N2MinHeapTiled) {
335 this->_minheap_faster_tiled_N2_cluster();
336 }
else if (_strategy == NlnN) {
337 this->_delaunay_cluster();
338 }
else if (_strategy == NlnNCam) {
339 this->_CP2DChan_cluster_2piMultD();
340 }
else if (_strategy == NlnN3pi || _strategy == NlnN4pi ) {
341 this->_delaunay_cluster();
342 }
else if (_strategy == N3Dumb ) {
343 this->_really_dumb_cluster();
344 }
else if (_strategy == N2PoorTiled) {
345 this->_tiled_N2_cluster();
346 }
else if (_strategy == NlnNCam4pi) {
347 this->_CP2DChan_cluster();
348 }
else if (_strategy == NlnNCam2pi2R) {
349 this->_CP2DChan_cluster_2pi2R();
352 err <<
"Unrecognised value for strategy: "<<_strategy;
353 throw Error(err.str());
360 bool ClusterSequence::_first_time =
true;
361 int ClusterSequence::_n_exclusive_warnings = 0;
367 return "FastJet version "+string(fastjet_version);
373 void ClusterSequence::print_banner() {
375 if (!_first_time) {
return;}
379 ostream * ostr = _fastjet_banner_ostr;
382 (*ostr) <<
"#--------------------------------------------------------------------------\n";
383 (*ostr) <<
"# FastJet release " << fastjet_version << endl;
384 (*ostr) <<
"# M. Cacciari, G.P. Salam and G. Soyez \n";
385 (*ostr) <<
"# A software package for jet finding and analysis at colliders \n";
386 (*ostr) <<
"# http://fastjet.fr \n";
388 (*ostr) <<
"# Please cite EPJC72(2012)1896 [arXiv:1111.6097] if you use this package\n";
389 (*ostr) <<
"# for scientific work and optionally PLB641(2006)57 [hep-ph/0512210]. \n";
391 (*ostr) <<
"# FastJet is provided without warranty under the terms of the GNU GPLv2.\n";
392 (*ostr) <<
"# It uses T. Chan's closest pair algorithm, S. Fortune's Voronoi code";
394 (*ostr) <<
",\n# CGAL ";
398 (*ostr) <<
"and 3rd party plugin jet algorithms. See COPYING file for details.\n";
399 (*ostr) <<
"#--------------------------------------------------------------------------\n";
407 const bool & writeout_combinations) {
409 _jet_def = jet_def_in;
410 _writeout_combinations = writeout_combinations;
414 _decant_options_partial();
419 void ClusterSequence::_decant_options_partial() {
423 _jet_algorithm = _jet_def.jet_algorithm();
424 _Rparam = _jet_def.R(); _R2 = _Rparam*_Rparam; _invR2 = 1.0/_R2;
425 _strategy = _jet_def.strategy();
428 _plugin_activated =
false;
432 _update_structure_use_count();
438 void ClusterSequence::_fill_initial_history () {
443 _jets.reserve(_jets.size()*2);
444 _history.reserve(_jets.size()*2);
448 for (
int i = 0; i < static_cast<int>(_jets.size()) ; i++) {
450 element.parent1 = InexistentParent;
451 element.
parent2 = InexistentParent;
452 element.
child = Invalid;
457 _history.push_back(element);
460 _jet_def.recombiner()->preprocess(_jets[i]);
463 _jets[i].set_cluster_hist_index(i);
464 _set_structure_shared_ptr(_jets[i]);
467 _Qtot += _jets[i].E();
469 _initial_n = _jets.size();
470 _deletes_self_when_unused =
false;
477 string ClusterSequence::strategy_string (
Strategy strategy_in)
const {
479 switch(strategy_in) {
481 strategy =
"NlnN";
break;
483 strategy =
"NlnN3pi";
break;
485 strategy =
"NlnN4pi";
break;
487 strategy =
"N2Plain";
break;
489 strategy =
"N2Tiled";
break;
491 strategy =
"N2MinHeapTiled";
break;
493 strategy =
"N2PoorTiled";
break;
495 strategy =
"N3Dumb";
break;
497 strategy =
"NlnNCam4pi";
break;
499 strategy =
"NlnNCam2pi2R";
break;
501 strategy =
"NlnNCam";
break;
503 strategy =
"plugin strategy";
break;
505 strategy =
"Unrecognized";
511 double ClusterSequence::jet_scale_for_algorithm(
516 double kt2=jet.
kt2();
517 return kt2 > 1e-300 ? 1.0/kt2 : 1e300;
519 double kt2 = jet.
kt2();
520 double p = jet_def().extra_param();
521 if (p <= 0 && kt2 < 1e-300) kt2 = 1e-300;
524 double kt2 = jet.
kt2();
525 double lim = _jet_def.extra_param();
526 if (kt2 < lim*lim && kt2 != 0.0) {
529 }
else {
throw Error(
"Unrecognised jet algorithm");}
590 if (will_delete_self_when_unused())
591 throw(
Error(
"cannot use CS::transfer_from_sequence after a call to delete_self_when_unused()"));
594 _jet_def = from_seq._jet_def ;
595 _writeout_combinations = from_seq._writeout_combinations ;
596 _initial_n = from_seq._initial_n ;
597 _Rparam = from_seq._Rparam ;
599 _invR2 = from_seq._invR2 ;
600 _strategy = from_seq._strategy ;
601 _jet_algorithm = from_seq._jet_algorithm ;
602 _plugin_activated = from_seq._plugin_activated ;
608 _jets = (*action_on_jets)(from_seq.
_jets);
610 _jets = from_seq.
_jets;
614 _extras = from_seq._extras;
617 if (_structure_shared_ptr()) {
621 if (_deletes_self_when_unused)
throw Error(
"transfer_from_sequence cannot be used for a cluster sequence that deletes self when unused");
632 _update_structure_use_count();
634 for (
unsigned int i=0; i<_jets.size(); i++){
637 _jets[i].set_cluster_hist_index(from_seq.
_jets[i].cluster_hist_index());
640 _set_structure_shared_ptr(_jets[i]);
648 void ClusterSequence::plugin_record_ij_recombination(
649 int jet_i,
int jet_j,
double dij,
650 const PseudoJet & newjet,
int & newjet_k) {
652 plugin_record_ij_recombination(jet_i, jet_j, dij, newjet_k);
655 int tmp_index = _jets[newjet_k].cluster_hist_index();
656 _jets[newjet_k] = newjet;
658 _set_structure_shared_ptr(_jets[newjet_k]);
664 vector<PseudoJet> ClusterSequence::inclusive_jets (
const double & ptmin)
const{
665 double dcut = ptmin*ptmin;
666 int i = _history.size() - 1;
667 vector<PseudoJet> jets_local;
673 if (_history[i].max_dij_so_far < dcut) {
break;}
674 if (_history[i].parent2 == BeamJet && _history[i].dij >= dcut) {
676 int parent1 = _history[i].parent1;
677 jets_local.push_back(_jets[_history[parent1].jetp_index]);}
685 if (_history[i].parent2 != BeamJet) {
break;}
686 int parent1 = _history[i].parent1;
687 const PseudoJet & jet = _jets[_history[parent1].jetp_index];
688 if (jet.
perp2() >= dcut) {jets_local.push_back(jet);}
701 if (_history[i].parent2 == BeamJet) {
702 int parent1 = _history[i].parent1;
703 const PseudoJet & jet = _jets[_history[parent1].jetp_index];
704 if (jet.
perp2() >= dcut) {jets_local.push_back(jet);}
708 }
else {
throw Error(
"cs::inclusive_jets(...): Unrecognized jet algorithm");}
716 int ClusterSequence::n_exclusive_jets (
const double & dcut)
const {
720 int i = _history.size() - 1;
722 if (_history[i].max_dij_so_far <= dcut) {
break;}
725 int stop_point = i + 1;
728 int njets = 2*_initial_n - stop_point;
735 vector<PseudoJet> ClusterSequence::exclusive_jets (
const double & dcut)
const {
736 int njets = n_exclusive_jets(dcut);
737 return exclusive_jets(njets);
744 vector<PseudoJet> ClusterSequence::exclusive_jets (
const int & njets)
const {
748 if (njets > _initial_n) {
750 err <<
"Requested " << njets <<
" exclusive jets, but there were only "
751 << _initial_n <<
" particles in the event";
752 throw Error(err.str());
755 return exclusive_jets_up_to(njets);
761 vector<PseudoJet> ClusterSequence::exclusive_jets_up_to (
const int & njets)
const {
773 (_jet_def.extra_param() <0)) &&
775 (!_jet_def.plugin()->exclusive_sequence_meaningful())) &&
776 (_n_exclusive_warnings < 5)) {
777 _n_exclusive_warnings++;
778 cerr <<
"FastJet WARNING: dcut and exclusive jets for jet-finders other than kt should be interpreted with care." << endl;
785 int stop_point = 2*_initial_n - njets;
787 if (stop_point < _initial_n) stop_point = _initial_n;
791 if (2*_initial_n != static_cast<int>(_history.size())) {
793 err <<
"2*_initial_n != _history.size() -- this endangers internal assumptions!\n";
794 throw Error(err.str());
803 vector<PseudoJet> jets_local;
804 for (
unsigned int i = stop_point; i < _history.size(); i++) {
805 int parent1 = _history[i].parent1;
806 if (parent1 < stop_point) {
807 jets_local.push_back(_jets[_history[parent1].jetp_index]);
809 int parent2 = _history[i].parent2;
810 if (parent2 < stop_point && parent2 > 0) {
811 jets_local.push_back(_jets[_history[parent2].jetp_index]);
817 if (
int(jets_local.size()) != min(_initial_n, njets)) {
819 err <<
"ClusterSequence::exclusive_jets: size of returned vector ("
820 <<jets_local.size()<<
") does not coincide with requested number of jets ("
822 throw Error(err.str());
831 double ClusterSequence::exclusive_dmerge (
const int & njets)
const {
833 if (njets >= _initial_n) {
return 0.0;}
834 return _history[2*_initial_n-njets-1].dij;
843 double ClusterSequence::exclusive_dmerge_max (
const int & njets)
const {
845 if (njets >= _initial_n) {
return 0.0;}
846 return _history[2*_initial_n-njets-1].max_dij_so_far;
854 std::vector<PseudoJet> ClusterSequence::exclusive_subjets
857 set<const history_element*> subhist;
861 get_subhist_set(subhist, jet, dcut, 0);
864 vector<PseudoJet> subjets;
865 subjets.reserve(subhist.size());
866 for (set<const history_element*>::iterator elem = subhist.begin();
867 elem != subhist.end(); elem++) {
868 subjets.push_back(_jets[(*elem)->jetp_index]);
877 int ClusterSequence::n_exclusive_subjets(
const PseudoJet & jet,
878 const double & dcut)
const {
879 set<const history_element*> subhist;
882 get_subhist_set(subhist, jet, dcut, 0);
883 return subhist.size();
890 std::vector<PseudoJet> ClusterSequence::exclusive_subjets
892 vector<PseudoJet> subjets = exclusive_subjets_up_to(jet, nsub);
893 if (
int(subjets.size()) < nsub) {
895 err <<
"Requested " << nsub <<
" exclusive subjets, but there were only "
896 << subjets.size() <<
" particles in the jet";
897 throw Error(err.str());
907 std::vector<PseudoJet> ClusterSequence::exclusive_subjets_up_to
910 set<const history_element*> subhist;
913 vector<PseudoJet> subjets;
914 if (nsub < 0)
throw Error(
"Requested a negative number of subjets. This is nonsensical.");
915 if (nsub == 0)
return subjets;
919 get_subhist_set(subhist, jet, -1.0, nsub);
922 subjets.reserve(subhist.size());
923 for (set<const history_element*>::iterator elem = subhist.begin();
924 elem != subhist.end(); elem++) {
925 subjets.push_back(_jets[(*elem)->jetp_index]);
936 double ClusterSequence::exclusive_subdmerge(
const PseudoJet & jet,
int nsub)
const {
937 set<const history_element*> subhist;
941 get_subhist_set(subhist, jet, -1.0, nsub);
943 set<const history_element*>::iterator highest = subhist.end();
947 return (*highest)->dij;
957 double ClusterSequence::exclusive_subdmerge_max(
const PseudoJet & jet,
int nsub)
const {
959 set<const history_element*> subhist;
963 get_subhist_set(subhist, jet, -1.0, nsub);
965 set<const history_element*>::iterator highest = subhist.end();
969 return (*highest)->max_dij_so_far;
981 void ClusterSequence::get_subhist_set(set<const history_element*> & subhist,
983 double dcut,
int maxjet)
const {
984 assert(contains(jet));
994 set<const history_element*>::iterator highest = subhist.end();
995 assert (highest != subhist.begin());
999 if (njet == maxjet)
break;
1001 if (elem->parent1 < 0)
break;
1006 subhist.erase(highest);
1007 subhist.insert(&(_history[elem->parent1]));
1008 subhist.insert(&(_history[elem->
parent2]));
1015 bool ClusterSequence::object_in_jet(
const PseudoJet &
object,
1020 assert(contains(
object) && contains(jet));
1022 const PseudoJet * this_object = &object;
1027 }
else if (has_child(*this_object, childp)) {
1028 this_object = childp;
1048 assert ((hist.parent1 >= 0 && hist.
parent2 >= 0) ||
1049 (hist.parent1 < 0 && hist.
parent2 < 0));
1051 if (hist.parent1 < 0) {
1056 parent1 = _jets[_history[hist.parent1].jetp_index];
1057 parent2 = _jets[_history[hist.
parent2].jetp_index];
1059 if (parent1.
perp2() < parent2.
perp2()) std::swap(parent1,parent2);
1079 bool res = has_child(jet, childp);
1096 if (hist.
child >= 0 && _history[hist.
child].jetp_index >= 0) {
1097 childp = &(_jets[_history[hist.
child].jetp_index]);
1117 if (hist.
child >= 0 && _history[hist.
child].parent2 >= 0) {
1122 partner = _jets[_history[child_hist.
parent2].jetp_index];
1125 partner = _jets[_history[child_hist.parent1].jetp_index];
1137 vector<PseudoJet> ClusterSequence::constituents (
const PseudoJet & jet)
const {
1138 vector<PseudoJet> subjets;
1139 add_constituents(jet, subjets);
1152 void ClusterSequence::print_jets_for_root(
const std::vector<PseudoJet> & jets_in,
1153 ostream & ostr)
const {
1154 for (
unsigned i = 0; i < jets_in.size(); i++) {
1156 << jets_in[i].px() <<
" "
1157 << jets_in[i].py() <<
" "
1158 << jets_in[i].pz() <<
" "
1159 << jets_in[i].E() << endl;
1160 vector<PseudoJet> cst = constituents(jets_in[i]);
1161 for (
unsigned j = 0; j < cst.size() ; j++) {
1162 ostr <<
" " << j <<
" "
1163 << cst[j].rap() <<
" "
1164 << cst[j].phi() <<
" "
1165 << cst[j].perp() << endl;
1167 ostr <<
"#END" << endl;
1171 void ClusterSequence::print_jets_for_root(
const std::vector<PseudoJet> & jets_in,
1172 const std::string & filename,
1173 const std::string & comment )
const {
1174 std::ofstream ostr(filename.c_str());
1175 if (comment !=
"") ostr <<
"# " << comment << endl;
1176 print_jets_for_root(jets_in, ostr);
1197 vector<int> ClusterSequence::particle_jet_indices(
1198 const vector<PseudoJet> & jets_in)
const {
1200 vector<int> indices(n_particles());
1203 for (
unsigned ipart = 0; ipart < n_particles(); ipart++)
1204 indices[ipart] = -1;
1208 for (
unsigned ijet = 0; ijet < jets_in.size(); ijet++) {
1210 vector<PseudoJet> jet_constituents(constituents(jets_in[ijet]));
1212 for (
unsigned ip = 0; ip < jet_constituents.size(); ip++) {
1216 unsigned iclust = jet_constituents[ip].cluster_hist_index();
1217 unsigned ipart = history()[iclust].jetp_index;
1218 indices[ipart] = ijet;
1228 void ClusterSequence::add_constituents (
1229 const PseudoJet & jet, vector<PseudoJet> & subjet_vector)
const {
1232 int parent1 = _history[i].parent1;
1233 int parent2 = _history[i].parent2;
1235 if (parent1 == InexistentParent) {
1241 subjet_vector.push_back(_jets[i]);
1246 add_constituents(_jets[_history[parent1].jetp_index], subjet_vector);
1249 if (parent2 != BeamJet) {
1250 add_constituents(_jets[_history[parent2].jetp_index], subjet_vector);
1258 void ClusterSequence::_add_step_to_history (
1259 const int & step_number,
const int & parent1,
1260 const int & parent2,
const int & jetp_index,
1261 const double & dij) {
1263 history_element element;
1264 element.parent1 = parent1;
1265 element.parent2 = parent2;
1266 element.jetp_index = jetp_index;
1267 element.child = Invalid;
1269 element.max_dij_so_far = max(dij,_history[_history.size()-1].max_dij_so_far);
1270 _history.push_back(element);
1272 int local_step = _history.size()-1;
1273 assert(local_step == step_number);
1275 assert(parent1 >= 0);
1276 _history[parent1].child = local_step;
1277 if (parent2 >= 0) {_history[parent2].child = local_step;}
1280 if (jetp_index != Invalid) {
1281 assert(jetp_index >= 0);
1283 _jets[jetp_index].set_cluster_hist_index(local_step);
1284 _set_structure_shared_ptr(_jets[jetp_index]);
1287 if (_writeout_combinations) {
1288 cout << local_step <<
": "
1289 << parent1 <<
" with " << parent2
1290 <<
"; y = "<< dij<<endl;
1303 vector<int> ClusterSequence::unique_history_order()
const {
1309 valarray<int> lowest_constituent(_history.size());
1310 int hist_n = _history.size();
1311 lowest_constituent = hist_n;
1312 for (
int i = 0; i < hist_n; i++) {
1314 lowest_constituent[i] = min(lowest_constituent[i],i);
1316 if (_history[i].child > 0) lowest_constituent[_history[i].child]
1317 = min(lowest_constituent[_history[i].child],lowest_constituent[i]);
1321 valarray<bool> extracted(_history.size()); extracted =
false;
1322 vector<int> unique_tree;
1323 unique_tree.reserve(_history.size());
1326 for (
unsigned i = 0; i < n_particles(); i++) {
1327 if (!extracted[i]) {
1328 unique_tree.push_back(i);
1329 extracted[i] =
true;
1330 _extract_tree_children(i, extracted, lowest_constituent, unique_tree);
1339 void ClusterSequence::_extract_tree_children(
1341 valarray<bool> & extracted,
1342 const valarray<int> & lowest_constituent,
1343 vector<int> & unique_tree)
const {
1344 if (!extracted[position]) {
1347 _extract_tree_parents(position,extracted,lowest_constituent,unique_tree);
1351 int child = _history[position].child;
1352 if (child >= 0) _extract_tree_children(child,extracted,lowest_constituent,unique_tree);
1358 vector<PseudoJet> ClusterSequence::unclustered_particles()
const {
1359 vector<PseudoJet> unclustered;
1360 for (
unsigned i = 0; i < n_particles() ; i++) {
1361 if (_history[i].child == Invalid)
1362 unclustered.push_back(_jets[_history[i].jetp_index]);
1372 vector<PseudoJet> ClusterSequence::childless_pseudojets()
const {
1373 vector<PseudoJet> unclustered;
1374 for (
unsigned i = 0; i < _history.size() ; i++) {
1375 if ((_history[i].child == Invalid) && (_history[i].parent2 != BeamJet))
1376 unclustered.push_back(_jets[_history[i].jetp_index]);
1387 bool ClusterSequence::contains(
const PseudoJet & jet)
const {
1398 void ClusterSequence::_extract_tree_parents(
1400 valarray<bool> & extracted,
1401 const valarray<int> & lowest_constituent,
1402 vector<int> & unique_tree)
const {
1404 if (!extracted[position]) {
1405 int parent1 = _history[position].parent1;
1406 int parent2 = _history[position].parent2;
1409 if (parent1 >= 0 && parent2 >= 0) {
1410 if (lowest_constituent[parent1] > lowest_constituent[parent2])
1411 std::swap(parent1, parent2);
1414 if (parent1 >= 0 && !extracted[parent1])
1415 _extract_tree_parents(parent1,extracted,lowest_constituent,unique_tree);
1416 if (parent2 >= 0 && !extracted[parent2])
1417 _extract_tree_parents(parent2,extracted,lowest_constituent,unique_tree);
1420 unique_tree.push_back(position);
1421 extracted[position] =
true;
1430 void ClusterSequence::_do_ij_recombination_step(
1431 const int & jet_i,
const int & jet_j,
1441 _jet_def.recombiner()->recombine(_jets[jet_i], _jets[jet_j], newjet);
1442 _jets.push_back(newjet);
1447 newjet_k = _jets.size()-1;
1450 int newstep_k = _history.size();
1452 _jets[newjet_k].set_cluster_hist_index(newstep_k);
1455 int hist_i = _jets[jet_i].cluster_hist_index();
1456 int hist_j = _jets[jet_j].cluster_hist_index();
1458 _add_step_to_history(newstep_k, min(hist_i, hist_j), max(hist_i,hist_j),
1467 void ClusterSequence::_do_iB_recombination_step(
1468 const int & jet_i,
const double & diB) {
1470 int newstep_k = _history.size();
1473 _add_step_to_history(newstep_k,_jets[jet_i].cluster_hist_index(),BeamJet,
1484 void ClusterSequence::_set_structure_shared_ptr(
PseudoJet & j) {
1488 _update_structure_use_count();
1493 void ClusterSequence::_update_structure_use_count() {
1496 _structure_use_count_after_construction = _structure_shared_ptr.use_count();
1503 void ClusterSequence::delete_self_when_unused() {
1511 int new_count = _structure_shared_ptr.use_count() - _structure_use_count_after_construction;
1512 if (new_count <= 0) {
1513 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");
1516 _structure_shared_ptr.set_count(new_count);
1517 _deletes_self_when_unused =
true;
1521 FASTJET_END_NAMESPACE
a version of cambridge with a special distance measure for particles whose pt is < extra_param() ...
like the k_t but with distance measures dij = min(kti^{2p},ktj^{2p}) Delta R_{ij}^2 / R^2 diB = 1/kti...
double kt2() const
returns the squared transverse momentum
Chan's closest pair method (in a variant with 4pi coverage), for use exclusively with the Cambridge a...
best of the NlnN variants – best overall for N>10^4.
legacy N ln N using 4pi coverage of cylinder
int cluster_hist_index() const
return the cluster_hist_index, intended to be used by clustering routines.
the e+e- genkt algorithm (R > 2 and p=1 gives ee_kt)
fastest from about 50..500
the longitudinally invariant kt algorithm
int jetp_index
index in _history where the current jet is recombined with another jet to form its child...
Contains any information related to the clustering that should be directly accessible to PseudoJet...
std::vector< PseudoJet > _jets
This contains the physical PseudoJets; for each PseudoJet one can find the corresponding position in ...
Chan's closest pair method (in a variant with 2pi+minimal extra variant), for use exclusively with th...
std::vector< history_element > _history
this vector will contain the branching history; for each stage, _history[i].jetp_index indicates wher...
class to provide facilities for giving warnings up to some maximum number of times and to provide glo...
any plugin algorithm supplied by the user
int parent2
index in _history where first parent of this jet was created (InexistentParent if this jet is an orig...
the plugin has been used...
virtual void set_associated_cs(const ClusterSequence *new_cs)
set the associated csw
double dij
index in the _jets vector where we will find the
void set_cluster_hist_index(const int index)
set the cluster_hist_index, intended to be used by clustering routines.
bool has_valid_cluster_sequence() const
returns true if this PseudoJet has an associated and still valid(ated) ClusterSequence.
worse even than the usual N^3 algorithms
base class corresponding to errors that can be thrown by FastJet
string fastjet_version_string()
return a string containing information about the release
int child
index in _history where second parent of this jet was created (InexistentParent if this jet is an ori...
automatic selection of the best (based on N)
a single element in the clustering history
like the k_t but with distance measures dij = min(1/kti^2,1/ktj^2) Delta R_{ij}^2 / R^2 diB = 1/kti^2...
const ClusterSequence * associated_cluster_sequence() const
get a (const) pointer to the parent ClusterSequence (NULL if inexistent)
void set_structure_shared_ptr(const SharedPtr< PseudoJetStructureBase > &structure)
set the associated structure
the longitudinally invariant variant of the cambridge algorithm (aka Aachen algoithm).
Strategy
the various options for the algorithmic strategy to adopt in clustering events with kt and cambridge ...
Class to contain pseudojets, including minimal information of use to jet-clustering routines...
double max_dij_so_far
the distance corresponding to the recombination at this stage of the clustering.
legacy N ln N using 3pi coverage of cylinder.
fastest form about 500..10^4
Chan's closest pair method (in a variant with 2pi+2R coverage), for use exclusively with the Cambridg...
double perp2() const
returns the squared transverse momentum
class that is intended to hold a full definition of the jet clusterer