30 #ifndef __FASTJET_CLUSTERSEQUENCE_HH__
31 #define __FASTJET_CLUSTERSEQUENCE_HH__
35 #include "fastjet/internal/DynamicNearestNeighbours.hh"
36 #include "fastjet/PseudoJet.hh"
43 #include "fastjet/Error.hh"
44 #include "fastjet/JetDefinition.hh"
45 #include "fastjet/SharedPtr.hh"
46 #include "fastjet/LimitedWarning.hh"
47 #include "fastjet/FunctionOfPseudoJet.hh"
48 #include "fastjet/ClusterSequenceStructure.hh"
50 FASTJET_BEGIN_NAMESPACE
54 class ClusterSequenceStructure;
86 const std::vector<L> & pseudojets,
88 const bool & writeout_combinations =
false);
92 transfer_from_sequence(cs);
106 std::vector<PseudoJet> inclusive_jets (
const double & ptmin = 0.0)
const;
111 int n_exclusive_jets (
const double & dcut)
const;
116 std::vector<PseudoJet> exclusive_jets (
const double & dcut)
const;
123 std::vector<PseudoJet> exclusive_jets (
const int & njets)
const;
130 std::vector<PseudoJet> exclusive_jets_up_to (
const int & njets)
const;
135 double exclusive_dmerge (
const int & njets)
const;
141 double exclusive_dmerge_max (
const int & njets)
const;
155 int njets = n_exclusive_jets_ycut(ycut);
156 return exclusive_jets(njets);
170 std::vector<PseudoJet> exclusive_subjets (
const PseudoJet & jet,
171 const double & dcut)
const;
176 int n_exclusive_subjets(
const PseudoJet & jet,
177 const double & dcut)
const;
184 std::vector<PseudoJet> exclusive_subjets (
const PseudoJet & jet,
192 std::vector<PseudoJet> exclusive_subjets_up_to (
const PseudoJet & jet,
199 double exclusive_subdmerge(
const PseudoJet & jet,
int nsub)
const;
206 double exclusive_subdmerge_max(
const PseudoJet & jet,
int nsub)
const;
213 double Q()
const {
return _Qtot;}
215 double Q2()
const {
return _Qtot*_Qtot;}
252 std::vector<PseudoJet> constituents (
const PseudoJet & jet)
const;
263 void print_jets_for_root(
const std::vector<PseudoJet> & jets,
264 std::ostream & ostr = std::cout)
const;
268 void print_jets_for_root(
const std::vector<PseudoJet> & jets,
269 const std::string & filename,
270 const std::string & comment =
"")
const;
277 void add_constituents (
const PseudoJet & jet,
278 std::vector<PseudoJet> & subjet_vector)
const;
287 std::string strategy_string (
Strategy strategy_in)
const;
305 void delete_self_when_unused();
312 void signal_imminent_self_deletion()
const;
317 double jet_scale_for_algorithm(
const PseudoJet & jet)
const;
329 void plugin_record_ij_recombination(
int jet_i,
int jet_j,
double dij,
331 assert(plugin_activated());
332 _do_ij_recombination_step(jet_i, jet_j, dij, newjet_k);
338 void plugin_record_ij_recombination(
int jet_i,
int jet_j,
double dij,
347 assert(plugin_activated());
348 _do_iB_recombination_step(jet_i, diB);
361 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";}
368 _extras.reset(extras_in.release());
385 assert(plugin_activated());
386 _simple_N2_cluster<GBJ>();
434 enum JetType {Invalid=-3, InexistentParent = -2, BeamJet = -1};
452 const std::vector<PseudoJet> & jets()
const;
465 const std::vector<history_element> & history()
const;
471 unsigned int n_particles()
const;
477 std::vector<int> particle_jet_indices(
const std::vector<PseudoJet> &)
const;
494 std::vector<int> unique_history_order()
const;
499 std::vector<PseudoJet> unclustered_particles()
const;
505 std::vector<PseudoJet> childless_pseudojets()
const;
512 bool contains(
const PseudoJet &
object)
const;
531 return _structure_shared_ptr;
544 static void print_banner();
559 static void set_fastjet_banner_stream(std::ostream * ostr) {_fastjet_banner_ostr = ostr;}
573 static std::ostream * _fastjet_banner_ostr;
583 template<
class L>
void _transfer_input_jets(
584 const std::vector<L> & pseudojets);
590 const bool & writeout_combinations);
596 void _initialise_and_run_no_decant();
608 const bool & writeout_combinations);
614 void _decant_options_partial();
619 void _fill_initial_history();
624 void _do_ij_recombination_step(
const int & jet_i,
const int & jet_j,
625 const double & dij,
int & newjet_k);
629 void _do_iB_recombination_step(
const int & jet_i,
const double & diB);
636 void _set_structure_shared_ptr(
PseudoJet & j);
640 void _update_structure_use_count();
662 void get_subhist_set(std::set<const history_element*> & subhist,
663 const PseudoJet & jet,
double dcut,
int maxjet)
const;
665 bool _writeout_combinations;
667 double _Rparam, _R2, _invR2;
673 int _structure_use_count_after_construction;
682 bool _plugin_activated;
686 void _really_dumb_cluster ();
687 void _delaunay_cluster ();
689 template<
class BJ>
void _simple_N2_cluster ();
690 void _tiled_N2_cluster ();
691 void _faster_tiled_N2_cluster ();
694 void _minheap_faster_tiled_N2_cluster();
698 void _CP2DChan_cluster();
699 void _CP2DChan_cluster_2pi2R ();
700 void _CP2DChan_cluster_2piMultD ();
701 void _CP2DChan_limited_cluster(
double D);
702 void _do_Cambridge_inclusive_jets();
705 void _fast_NsqrtN_cluster();
707 void _add_step_to_history(
const int & step_number,
const int & parent1,
708 const int & parent2,
const int & jetp_index,
713 void _extract_tree_children(
int pos, std::valarray<bool> &,
714 const std::valarray<int> &, std::vector<int> &)
const;
718 void _extract_tree_parents (
int pos, std::valarray<bool> &,
719 const std::valarray<int> &, std::vector<int> &)
const;
723 typedef std::pair<int,int> TwoVertices;
724 typedef std::pair<double,TwoVertices> DijEntry;
725 typedef std::multimap<double,TwoVertices> DistMap;
728 void _add_ktdistance_to_map(
const int & ii,
730 const DynamicNearestNeighbours * DNN);
734 static bool _first_time;
739 static int _n_exclusive_warnings;
752 double eta, phi, kt2, NN_dist;
762 double eta, phi, kt2, NN_dist;
763 TiledJet * NN, *previous, * next;
764 int _jets_index, tile_index, diJ_posn;
769 inline void label_minheap_update_needed() {diJ_posn = 1;}
770 inline void label_minheap_update_done() {diJ_posn = 0;}
771 inline bool minheap_update_needed()
const {
return diJ_posn==1;}
779 template <
class J>
void _bj_set_jetinfo( J *
const jet,
780 const int _jets_index)
const;
784 void _bj_remove_from_tiles( TiledJet *
const jet)
const;
787 template <
class J>
double _bj_dist(
const J *
const jeta,
788 const J *
const jetb)
const;
792 template <
class J>
double _bj_diJ(
const J *
const jeta)
const;
797 template <
class J>
inline J * _bj_of_hindex(
798 const int hist_index,
799 J *
const head, J *
const tail)
802 for(res = head; res<tail; res++) {
803 if (_jets[res->_jets_index].cluster_hist_index() == hist_index) {
break;}
816 template <
class J>
void _bj_set_NN_nocross(J *
const jeta,
817 J *
const head,
const J *
const tail)
const;
823 template <
class J>
void _bj_set_NN_crosscheck(J *
const jeta,
824 J *
const head,
const J *
const tail)
const;
830 static const int n_tile_neighbours = 9;
836 Tile * begin_tiles[n_tile_neighbours];
838 Tile ** surrounding_tiles;
848 std::vector<Tile> _tiles;
849 double _tiles_eta_min, _tiles_eta_max;
850 double _tile_size_eta, _tile_size_phi;
851 int _n_tiles_phi,_tiles_ieta_min,_tiles_ieta_max;
855 inline int _tile_index (
int ieta,
int iphi)
const {
858 return (ieta-_tiles_ieta_min)*_n_tiles_phi
859 + (iphi+_n_tiles_phi) % _n_tiles_phi;
864 int _tile_index(
const double & eta,
const double & phi)
const;
865 void _tj_set_jetinfo ( TiledJet *
const jet,
const int _jets_index);
866 void _bj_remove_from_tiles(TiledJet *
const jet);
867 void _initialise_tiles();
868 void _print_tiles(TiledJet * briefjets )
const;
869 void _add_neighbours_to_tile_union(
const int tile_index,
870 std::vector<int> & tile_union,
int & n_near_tiles)
const;
871 void _add_untagged_neighbours_to_tile_union(
const int tile_index,
872 std::vector<int> & tile_union,
int & n_near_tiles);
891 void _simple_N2_cluster_BriefJet();
893 void _simple_N2_cluster_EEBriefJet();
904 template<
class L>
void ClusterSequence::_transfer_input_jets(
905 const std::vector<L> & pseudojets) {
909 _jets.reserve(pseudojets.size()*2);
915 for (
unsigned int i = 0; i < pseudojets.size(); i++) {
916 _jets.push_back(pseudojets[i]);}
940 template<
class L> ClusterSequence::ClusterSequence (
941 const std::vector<L> & pseudojets,
943 const bool & writeout_combinations) :
944 _jet_def(jet_def_in), _writeout_combinations(writeout_combinations),
955 _initialise_and_run_no_decant();
972 template <
class J>
inline void ClusterSequence::_bj_set_jetinfo(
973 J *
const jetA,
const int _jets_index)
const {
974 jetA->eta =
_jets[_jets_index].rap();
975 jetA->phi =
_jets[_jets_index].phi_02pi();
977 jetA->_jets_index = _jets_index;
987 template <
class J>
inline double ClusterSequence::_bj_dist(
988 const J *
const jetA,
const J *
const jetB)
const {
989 double dphi = std::abs(jetA->phi - jetB->phi);
990 double deta = (jetA->eta - jetB->eta);
991 if (dphi > pi) {dphi = twopi - dphi;}
992 return dphi*dphi + deta*deta;
996 template <
class J>
inline double ClusterSequence::_bj_diJ(
const J *
const jet)
const {
997 double kt2 = jet->kt2;
998 if (jet->NN != NULL) {
if (jet->NN->kt2 < kt2) {kt2 = jet->NN->kt2;}}
999 return jet->NN_dist * kt2;
1006 template <
class J>
inline void ClusterSequence::_bj_set_NN_nocross(
1007 J *
const jet, J *
const head,
const J *
const tail)
const {
1008 double NN_dist = _R2;
1011 for (J * jetB = head; jetB != jet; jetB++) {
1012 double dist = _bj_dist(jet,jetB);
1013 if (dist < NN_dist) {
1020 for (J * jetB = jet+1; jetB != tail; jetB++) {
1021 double dist = _bj_dist(jet,jetB);
1022 if (dist < NN_dist) {
1029 jet->NN_dist = NN_dist;
1034 template <
class J>
inline void ClusterSequence::_bj_set_NN_crosscheck(J *
const jet,
1035 J *
const head,
const J *
const tail)
const {
1036 double NN_dist = _R2;
1038 for (J * jetB = head; jetB != tail; jetB++) {
1039 double dist = _bj_dist(jet,jetB);
1040 if (dist < NN_dist) {
1044 if (dist < jetB->NN_dist) {
1045 jetB->NN_dist = dist;
1050 jet->NN_dist = NN_dist;
1053 FASTJET_END_NAMESPACE
1055 #endif // __FASTJET_CLUSTERSEQUENCE_HH__