1 #ifndef __FASTJET_CLUSTERSEQUENCE_HH__
2 #define __FASTJET_CLUSTERSEQUENCE_HH__
37 #include "fastjet/PseudoJet.hh"
44 #include "fastjet/Error.hh"
45 #include "fastjet/JetDefinition.hh"
46 #include "fastjet/SharedPtr.hh"
47 #include "fastjet/LimitedWarning.hh"
48 #include "fastjet/FunctionOfPseudoJet.hh"
49 #include "fastjet/ClusterSequenceStructure.hh"
50 #include "fastjet/internal/thread_safety_helpers.hh"
51 #include "fastjet/internal/deprecated.hh"
53 FASTJET_BEGIN_NAMESPACE
57 class ClusterSequenceStructure;
58 class DynamicNearestNeighbours;
75 const std::vector<L> & pseudojets,
77 const bool & writeout_combinations =
false);
81 transfer_from_sequence(cs);
98 std::vector<PseudoJet> inclusive_jets (
const double ptmin = 0.0)
const;
103 int n_exclusive_jets (
const double dcut)
const;
108 std::vector<PseudoJet> exclusive_jets (
const double dcut)
const;
115 std::vector<PseudoJet> exclusive_jets (
const int njets)
const;
122 std::vector<PseudoJet> exclusive_jets_up_to (
const int njets)
const;
127 double exclusive_dmerge (
const int njets)
const;
133 double exclusive_dmerge_max (
const int njets)
const;
147 int njets = n_exclusive_jets_ycut(ycut);
148 return exclusive_jets(njets);
162 std::vector<PseudoJet> exclusive_subjets (
const PseudoJet & jet,
163 const double dcut)
const;
168 int n_exclusive_subjets(
const PseudoJet & jet,
169 const double dcut)
const;
176 std::vector<PseudoJet> exclusive_subjets (
const PseudoJet & jet,
184 std::vector<PseudoJet> exclusive_subjets_up_to (
const PseudoJet & jet,
191 double exclusive_subdmerge(
const PseudoJet & jet,
int nsub)
const;
198 double exclusive_subdmerge_max(
const PseudoJet & jet,
int nsub)
const;
205 double Q()
const {
return _Qtot;}
207 double Q2()
const {
return _Qtot*_Qtot;}
244 std::vector<PseudoJet> constituents (
const PseudoJet & jet)
const;
256 std::ostream & ostr = std::cout)
const;
260 void print_jets_for_root(
const std::vector<PseudoJet> & jets,
261 const std::string & filename,
262 const std::string & comment =
"")
const;
269 void add_constituents (
const PseudoJet & jet,
270 std::vector<PseudoJet> & subjet_vector)
const;
279 std::string strategy_string (
Strategy strategy_in)
const;
296 void delete_self_when_unused();
307 void signal_imminent_self_deletion()
const;
314 double jet_scale_for_algorithm(
const PseudoJet & jet)
const;
328 assert(plugin_activated());
329 _do_ij_recombination_step(jet_i, jet_j, dij, newjet_k);
335 void plugin_record_ij_recombination(
int jet_i,
int jet_j,
double dij,
344 assert(plugin_activated());
345 _do_iB_recombination_step(jet_i, diB);
358 virtual std::string description()
const {
return "This is a dummy extras class that contains no extra information! Derive from it if you want to use it to provide extra information from a plugin jet finder";}
366 _extras.reset(extras_in);
375 #ifdef FASTJET_HAVE_AUTO_PTR_INTERFACE
376 FASTJET_DEPRECATED_MSG(
"Please use ClusterSequence::plugin_associate_extras(Extras * extras_in)) instead",
378 _extras.reset(extras_in.release());
397 assert(plugin_activated());
398 _simple_N2_cluster<GBJ>();
446 enum JetType {Invalid=-3, InexistentParent = -2, BeamJet = -1};
464 const std::vector<PseudoJet> & jets()
const;
477 const std::vector<history_element> & history()
const;
483 unsigned int n_particles()
const;
489 std::vector<int> particle_jet_indices(
const std::vector<PseudoJet> &)
const;
506 std::vector<int> unique_history_order()
const;
511 std::vector<PseudoJet> unclustered_particles()
const;
517 std::vector<PseudoJet> childless_pseudojets()
const;
524 bool contains(
const PseudoJet &
object)
const;
543 return _structure_shared_ptr;
556 static void print_banner();
571 static void set_fastjet_banner_stream(std::ostream * ostr) {_fastjet_banner_ostr = ostr;}
585 #ifdef FASTJET_HAVE_LIMITED_THREAD_SAFETY
586 static std::atomic<std::ostream*> _fastjet_banner_ostr;
588 static std::ostream * _fastjet_banner_ostr;
589 #endif // FASTJET_HAVE_LIMITED_THREAD_SAFETY
598 template<
class L>
void _transfer_input_jets(
599 const std::vector<L> & pseudojets);
605 const bool & writeout_combinations);
611 void _initialise_and_run_no_decant();
623 const bool & writeout_combinations);
629 void _decant_options_partial();
634 void _fill_initial_history();
639 void _do_ij_recombination_step(
const int jet_i,
const int jet_j,
640 const double dij,
int & newjet_k);
644 void _do_iB_recombination_step(
const int jet_i,
const double diB);
651 void _set_structure_shared_ptr(
PseudoJet & j);
655 void _update_structure_use_count();
669 _Parabola(
double a,
double b,
double c) : _a(a), _b(b), _c(c) {}
670 inline double operator()(
const double R)
const {
return _c*(_a*R*R + _b*R + 1);}
681 _Line(
double a,
double b) : _a(a), _b(b) {}
682 inline double operator()(
const double R)
const {
return _a*R + _b;}
706 void get_subhist_set(std::set<const history_element*> & subhist,
707 const PseudoJet & jet,
double dcut,
int maxjet)
const;
709 bool _writeout_combinations;
711 double _Rparam, _R2, _invR2;
717 int _structure_use_count_after_construction;
722 #ifdef FASTJET_HAVE_THREAD_SAFETY
723 mutable std::atomic<bool> _deletes_self_when_unused;
730 bool _plugin_activated;
733 void _really_dumb_cluster ();
734 void _delaunay_cluster ();
736 template<
class BJ>
void _simple_N2_cluster ();
737 void _tiled_N2_cluster ();
738 void _faster_tiled_N2_cluster ();
741 void _minheap_faster_tiled_N2_cluster();
745 void _CP2DChan_cluster();
746 void _CP2DChan_cluster_2pi2R ();
747 void _CP2DChan_cluster_2piMultD ();
748 void _CP2DChan_limited_cluster(
double D);
749 void _do_Cambridge_inclusive_jets();
752 void _fast_NsqrtN_cluster();
754 void _add_step_to_history(
756 const int parent2,
const int jetp_index,
761 void _extract_tree_children(
int pos, std::valarray<bool> &,
762 const std::valarray<int> &, std::vector<int> &)
const;
766 void _extract_tree_parents (
int pos, std::valarray<bool> &,
767 const std::valarray<int> &, std::vector<int> &)
const;
771 typedef std::pair<int,int> TwoVertices;
772 typedef std::pair<double,TwoVertices> DijEntry;
773 typedef std::multimap<double,TwoVertices> DistMap;
776 void _add_ktdistance_to_map(
const int ii,
778 const DynamicNearestNeighbours * DNN);
798 double eta, phi, kt2, NN_dist;
807 double eta, phi, kt2, NN_dist;
809 int _jets_index, tile_index, diJ_posn;
814 inline void label_minheap_update_needed() {diJ_posn = 1;}
815 inline void label_minheap_update_done() {diJ_posn = 0;}
816 inline bool minheap_update_needed()
const {
return diJ_posn==1;}
824 template <
class J>
void _bj_set_jetinfo( J *
const jet,
825 const int _jets_index)
const;
829 void _bj_remove_from_tiles( TiledJet *
const jet)
const;
832 template <
class J>
double _bj_dist(
const J *
const jeta,
833 const J *
const jetb)
const;
837 template <
class J>
double _bj_diJ(
const J *
const jeta)
const;
842 template <
class J>
inline J * _bj_of_hindex(
843 const int hist_index,
844 J *
const head, J *
const tail)
847 for(res = head; res<tail; res++) {
848 if (_jets[res->_jets_index].cluster_hist_index() == hist_index) {
break;}
861 template <
class J>
void _bj_set_NN_nocross(J *
const jeta,
862 J *
const head,
const J *
const tail)
const;
868 template <
class J>
void _bj_set_NN_crosscheck(J *
const jeta,
869 J *
const head,
const J *
const tail)
const;
875 static const int n_tile_neighbours = 9;
881 Tile * begin_tiles[n_tile_neighbours];
883 Tile ** surrounding_tiles;
893 std::vector<Tile> _tiles;
894 double _tiles_eta_min, _tiles_eta_max;
895 double _tile_size_eta, _tile_size_phi;
896 int _n_tiles_phi,_tiles_ieta_min,_tiles_ieta_max;
900 inline int _tile_index (
int ieta,
int iphi)
const {
903 return (ieta-_tiles_ieta_min)*_n_tiles_phi
904 + (iphi+_n_tiles_phi) % _n_tiles_phi;
909 int _tile_index(
const double eta,
const double phi)
const;
910 void _tj_set_jetinfo ( TiledJet *
const jet,
const int _jets_index);
911 void _bj_remove_from_tiles(TiledJet *
const jet);
912 void _initialise_tiles();
913 void _print_tiles(TiledJet * briefjets )
const;
914 void _add_neighbours_to_tile_union(
const int tile_index,
915 std::vector<int> & tile_union,
int & n_near_tiles)
const;
916 void _add_untagged_neighbours_to_tile_union(
const int tile_index,
917 std::vector<int> & tile_union,
int & n_near_tiles);
935 void _simple_N2_cluster_BriefJet();
937 void _simple_N2_cluster_EEBriefJet();
948 template<
class L>
void ClusterSequence::_transfer_input_jets(
949 const std::vector<L> & pseudojets) {
953 _jets.reserve(pseudojets.size()*2);
959 for (
unsigned int i = 0; i < pseudojets.size(); i++) {
960 _jets.push_back(pseudojets[i]);}
984 template<
class L> ClusterSequence::ClusterSequence (
985 const std::vector<L> & pseudojets,
987 const bool & writeout_combinations) :
988 _jet_def(jet_def_in), _writeout_combinations(writeout_combinations),
999 _initialise_and_run_no_decant();
1024 std::vector<PseudoJet> jets;
1033 if (jets.size() != 0) {
1045 template <
class J>
inline void ClusterSequence::_bj_set_jetinfo(
1046 J *
const jetA,
const int _jets_index)
const {
1047 jetA->eta =
_jets[_jets_index].rap();
1048 jetA->phi =
_jets[_jets_index].phi_02pi();
1050 jetA->_jets_index = _jets_index;
1052 jetA->NN_dist = _R2;
1060 template <
class J>
inline double ClusterSequence::_bj_dist(
1061 const J *
const jetA,
const J *
const jetB)
const {
1063 #ifndef FASTJET_NEW_DELTA_PHI
1065 double dphi = std::abs(jetA->phi - jetB->phi);
1066 double deta = (jetA->eta - jetB->eta);
1067 if (dphi > pi) {dphi = twopi - dphi;}
1070 double dphi = pi-std::abs(pi-std::abs(jetA->phi - jetB->phi));
1071 double deta = (jetA->eta - jetB->eta);
1073 return dphi*dphi + deta*deta;
1077 template <
class J>
inline double ClusterSequence::_bj_diJ(
const J *
const jet)
const {
1078 double kt2 = jet->kt2;
1079 if (jet->NN != NULL) {
if (jet->NN->kt2 < kt2) {kt2 = jet->NN->kt2;}}
1080 return jet->NN_dist * kt2;
1087 template <
class J>
inline void ClusterSequence::_bj_set_NN_nocross(
1088 J *
const jet, J *
const head,
const J *
const tail)
const {
1089 double NN_dist = _R2;
1092 for (J * jetB = head; jetB != jet; jetB++) {
1093 double dist = _bj_dist(jet,jetB);
1094 if (dist < NN_dist) {
1101 for (J * jetB = jet+1; jetB != tail; jetB++) {
1102 double dist = _bj_dist(jet,jetB);
1103 if (dist < NN_dist) {
1110 jet->NN_dist = NN_dist;
1115 template <
class J>
inline void ClusterSequence::_bj_set_NN_crosscheck(J *
const jet,
1116 J *
const head,
const J *
const tail)
const {
1117 double NN_dist = _R2;
1119 for (J * jetB = head; jetB != tail; jetB++) {
1120 double dist = _bj_dist(jet,jetB);
1121 if (dist < NN_dist) {
1125 if (dist < jetB->NN_dist) {
1126 jetB->NN_dist = dist;
1131 jet->NN_dist = NN_dist;
1134 FASTJET_END_NAMESPACE
1136 #endif // __FASTJET_CLUSTERSEQUENCE_HH__