00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013
00014
00015
00016
00017
00018
00019
00020
00021
00022
00023
00024
00025
00026
00027
00028
00029
00030
00031
00032
00033
00034
00035
00044
00045
00046 #ifndef __FASTJET_CLUSTERSEQUENCE_HH__
00047 #define __FASTJET_CLUSTERSEQUENCE_HH__
00048
00049 #include<vector>
00050 #include<map>
00051 #include "fastjet/internal/DynamicNearestNeighbours.hh"
00052 #include "fastjet/PseudoJet.hh"
00053 #include<memory>
00054 #include<cassert>
00055 #include<iostream>
00056 #include<string>
00057 #include<set>
00058 #include<cmath>
00059 #include "fastjet/Error.hh"
00060 #include "fastjet/JetDefinition.hh"
00061
00062 FASTJET_BEGIN_NAMESPACE
00063
00064
00066 class ClusterSequence {
00067
00068
00069 public:
00070
00072 ClusterSequence () {}
00073
00083 template<class L> ClusterSequence (const std::vector<L> & pseudojets,
00084 const double & R = 1.0,
00085 const Strategy & strategy = Best,
00086 const bool & writeout_combinations = false);
00087
00088
00092 template<class L> ClusterSequence (
00093 const std::vector<L> & pseudojets,
00094 const JetDefinition & jet_def,
00095 const bool & writeout_combinations = false);
00096
00097
00098
00099 virtual ~ClusterSequence ();
00100
00101
00102
00103
00104
00108 std::vector<PseudoJet> inclusive_jets (const double & ptmin = 0.0) const;
00109
00113 int n_exclusive_jets (const double & dcut) const;
00114
00118 std::vector<PseudoJet> exclusive_jets (const double & dcut) const;
00119
00122 std::vector<PseudoJet> exclusive_jets (const int & njets) const;
00123
00126 double exclusive_dmerge (const int & njets) const;
00127
00132 double exclusive_dmerge_max (const int & njets) const;
00133
00136 double exclusive_ymerge (int njets) const {return exclusive_dmerge(njets) / Q2();}
00137
00139 double exclusive_ymerge_max (int njets) const {return exclusive_dmerge_max(njets)/Q2();}
00140
00142 int n_exclusive_jets_ycut (double ycut) const {return n_exclusive_jets(ycut*Q2());}
00143
00145 std::vector<PseudoJet> exclusive_jets_ycut (double ycut) const {
00146 int njets = n_exclusive_jets_ycut(ycut);
00147 return exclusive_jets(njets);
00148 }
00149
00150
00151
00152
00161 std::vector<PseudoJet> exclusive_subjets (const PseudoJet & jet,
00162 const double & dcut) const;
00163
00167 int n_exclusive_subjets(const PseudoJet & jet,
00168 const double & dcut) const;
00169
00175 std::vector<PseudoJet> exclusive_subjets (const PseudoJet & jet,
00176 int nsub) const;
00177
00180 double exclusive_subdmerge(const PseudoJet & jet, int nsub) const;
00181
00185 double exclusive_subdmerge_max(const PseudoJet & jet, int nsub) const;
00186
00187
00188
00189
00190
00192 double Q() const {return _Qtot;}
00194 double Q2() const {return _Qtot*_Qtot;}
00195
00206 bool object_in_jet(const PseudoJet & object, const PseudoJet & jet) const;
00207
00213 bool has_parents(const PseudoJet & jet, PseudoJet & parent1,
00214 PseudoJet & parent2) const;
00215
00218 bool has_child(const PseudoJet & jet, PseudoJet & child) const;
00219
00222 bool has_child(const PseudoJet & jet, const PseudoJet * & childp) const;
00223
00227 bool has_partner(const PseudoJet & jet, PseudoJet & partner) const;
00228
00229
00231 std::vector<PseudoJet> constituents (const PseudoJet & jet) const;
00232
00233
00242 void print_jets_for_root(const std::vector<PseudoJet> & jets,
00243 std::ostream & ostr = std::cout) const;
00244
00247 void print_jets_for_root(const std::vector<PseudoJet> & jets,
00248 const std::string & filename,
00249 const std::string & comment = "") const;
00250
00251
00252
00253
00254
00256 void add_constituents (const PseudoJet & jet,
00257 std::vector<PseudoJet> & subjet_vector) const;
00258
00260 inline Strategy strategy_used () const {return _strategy;}
00261 std::string strategy_string () const;
00262
00264 const JetDefinition & jet_def() const {return _jet_def;}
00265
00269 double jet_scale_for_algorithm(const PseudoJet & jet) const;
00270
00271
00272
00273
00279 void plugin_record_ij_recombination(int jet_i, int jet_j, double dij,
00280 int & newjet_k) {
00281 assert(plugin_activated());
00282 _do_ij_recombination_step(jet_i, jet_j, dij, newjet_k);
00283 }
00284
00288 void plugin_record_ij_recombination(int jet_i, int jet_j, double dij,
00289 const PseudoJet & newjet,
00290 int & newjet_k);
00291
00296 void plugin_record_iB_recombination(int jet_i, double diB) {
00297 assert(plugin_activated());
00298 _do_iB_recombination_step(jet_i, diB);
00299 }
00300
00304 class Extras {
00305 public:
00306 virtual ~Extras() {}
00307 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";}
00308 };
00309
00312 inline void plugin_associate_extras(std::auto_ptr<Extras> extras_in) {
00313 _extras = extras_in;
00314 }
00315
00317 inline bool plugin_activated() const {return _plugin_activated;}
00318
00320 const Extras * extras() const {return _extras.get();}
00321
00329 template<class GBJ> void plugin_simple_N2_cluster () {
00330 assert(plugin_activated());
00331 _simple_N2_cluster<GBJ>();
00332 }
00333
00334
00335
00336
00337
00338
00339
00340
00341
00342
00343
00344
00345
00346
00347
00348
00349
00350
00351
00352
00353
00354
00355
00356
00357
00358
00359
00360
00361
00362
00363
00364
00365
00366
00367 public:
00371 static void set_jet_algorithm (JetAlgorithm jet_algorithm) {_default_jet_algorithm = jet_algorithm;}
00373 static void set_jet_finder (JetAlgorithm jet_algorithm) {_default_jet_algorithm = jet_algorithm;}
00374
00375
00378 struct history_element{
00379 int parent1;
00380
00381
00382
00383 int parent2;
00384
00385
00386
00387
00388
00389 int child;
00390
00391
00392
00393
00394 int jetp_index;
00395
00396
00397
00398
00399
00400
00401 double dij;
00402
00403
00404 double max_dij_so_far;
00405
00406 };
00407
00408 enum JetType {Invalid=-3, InexistentParent = -2, BeamJet = -1};
00409
00426 const std::vector<PseudoJet> & jets() const;
00427
00439 const std::vector<history_element> & history() const;
00440
00445 unsigned int n_particles() const;
00446
00451 std::vector<int> particle_jet_indices(const std::vector<PseudoJet> &) const;
00452
00468 std::vector<int> unique_history_order() const;
00469
00473 std::vector<PseudoJet> unclustered_particles() const;
00474
00478 void transfer_from_sequence(ClusterSequence & from_seq);
00479
00480
00481 protected:
00482 static JetAlgorithm _default_jet_algorithm;
00483 JetDefinition _jet_def;
00484
00487 bool _potentially_valid(const PseudoJet & jet) const {
00488 return jet.cluster_hist_index() >= 0
00489 && jet.cluster_hist_index() < int(_history.size());
00490 }
00491
00494 template<class L> void _transfer_input_jets(
00495 const std::vector<L> & pseudojets);
00496
00500 void _initialise_and_run (const JetDefinition & jet_def,
00501 const bool & writeout_combinations);
00502
00506 void _initialise_and_run (const double & R,
00507 const Strategy & strategy,
00508 const bool & writeout_combinations);
00509
00512 void _decant_options(const JetDefinition & jet_def,
00513 const bool & writeout_combinations);
00514
00518 void _fill_initial_history();
00519
00523 void _do_ij_recombination_step(const int & jet_i, const int & jet_j,
00524 const double & dij, int & newjet_k);
00525
00528 void _do_iB_recombination_step(const int & jet_i, const double & diB);
00529
00530
00534 std::vector<PseudoJet> _jets;
00535
00536
00540 std::vector<history_element> _history;
00541
00550 void get_subhist_set(std::set<const history_element*> & subhist,
00551 const PseudoJet & jet, double dcut, int maxjet) const;
00552
00553 bool _writeout_combinations;
00554 int _initial_n;
00555 double _Rparam, _R2, _invR2;
00556 double _Qtot;
00557 Strategy _strategy;
00558 JetAlgorithm _jet_algorithm;
00559
00560
00561 private:
00562
00563 bool _plugin_activated;
00564 std::auto_ptr<Extras> _extras;
00565
00566 void _really_dumb_cluster ();
00567 void _delaunay_cluster ();
00568
00569 template<class BJ> void _simple_N2_cluster ();
00570 void _tiled_N2_cluster ();
00571 void _faster_tiled_N2_cluster ();
00572
00573
00574 void _minheap_faster_tiled_N2_cluster();
00575
00576
00577
00578 void _CP2DChan_cluster();
00579 void _CP2DChan_cluster_2pi2R ();
00580 void _CP2DChan_cluster_2piMultD ();
00581 void _CP2DChan_limited_cluster(double D);
00582 void _do_Cambridge_inclusive_jets();
00583
00584 void _add_step_to_history(const int & step_number, const int & parent1,
00585 const int & parent2, const int & jetp_index,
00586 const double & dij);
00587
00590 void _extract_tree_children(int pos, std::valarray<bool> &,
00591 const std::valarray<int> &, std::vector<int> &) const;
00592
00595 void _extract_tree_parents (int pos, std::valarray<bool> &,
00596 const std::valarray<int> &, std::vector<int> &) const;
00597
00598
00599
00600 typedef std::pair<int,int> TwoVertices;
00601 typedef std::pair<double,TwoVertices> DijEntry;
00602 typedef std::multimap<double,TwoVertices> DistMap;
00603
00605 void _add_ktdistance_to_map(const int & ii,
00606 DistMap & DijMap,
00607 const DynamicNearestNeighbours * DNN);
00608
00610 void _print_banner();
00612 static bool _first_time;
00613
00617 static int _n_exclusive_warnings;
00618
00619
00624 struct BriefJet {
00625 double eta, phi, kt2, NN_dist;
00626 BriefJet * NN;
00627 int _jets_index;
00628 };
00629
00630
00633 class TiledJet {
00634 public:
00635 double eta, phi, kt2, NN_dist;
00636 TiledJet * NN, *previous, * next;
00637 int _jets_index, tile_index, diJ_posn;
00638
00639
00640
00641
00642 inline void label_minheap_update_needed() {diJ_posn = 1;}
00643 inline void label_minheap_update_done() {diJ_posn = 0;}
00644 inline bool minheap_update_needed() const {return diJ_posn==1;}
00645 };
00646
00647
00648
00649
00652 template <class J> void _bj_set_jetinfo( J * const jet,
00653 const int _jets_index) const;
00654
00657 void _bj_remove_from_tiles( TiledJet * const jet) const;
00658
00660 template <class J> double _bj_dist(const J * const jeta,
00661 const J * const jetb) const;
00662
00663
00664
00665 template <class J> double _bj_diJ(const J * const jeta) const;
00666
00670 template <class J> inline J * _bj_of_hindex(
00671 const int hist_index,
00672 J * const head, J * const tail)
00673 const {
00674 J * res;
00675 for(res = head; res<tail; res++) {
00676 if (_jets[res->_jets_index].cluster_hist_index() == hist_index) {break;}
00677 }
00678 return res;
00679 }
00680
00681
00682
00683
00684
00689 template <class J> void _bj_set_NN_nocross(J * const jeta,
00690 J * const head, const J * const tail) const;
00691
00696 template <class J> void _bj_set_NN_crosscheck(J * const jeta,
00697 J * const head, const J * const tail) const;
00698
00699
00700
00703 static const int n_tile_neighbours = 9;
00704
00707 struct Tile {
00709 Tile * begin_tiles[n_tile_neighbours];
00711 Tile ** surrounding_tiles;
00713 Tile ** RH_tiles;
00715 Tile ** end_tiles;
00717 TiledJet * head;
00719 bool tagged;
00720 };
00721 std::vector<Tile> _tiles;
00722 double _tiles_eta_min, _tiles_eta_max;
00723 double _tile_size_eta, _tile_size_phi;
00724 int _n_tiles_phi,_tiles_ieta_min,_tiles_ieta_max;
00725
00726
00727
00728 inline int _tile_index (int ieta, int iphi) const {
00729
00730
00731 return (ieta-_tiles_ieta_min)*_n_tiles_phi
00732 + (iphi+_n_tiles_phi) % _n_tiles_phi;
00733 }
00734
00735
00736
00737 int _tile_index(const double & eta, const double & phi) const;
00738 void _tj_set_jetinfo ( TiledJet * const jet, const int _jets_index);
00739 void _bj_remove_from_tiles(TiledJet * const jet);
00740 void _initialise_tiles();
00741 void _print_tiles(TiledJet * briefjets ) const;
00742 void _add_neighbours_to_tile_union(const int tile_index,
00743 std::vector<int> & tile_union, int & n_near_tiles) const;
00744 void _add_untagged_neighbours_to_tile_union(const int tile_index,
00745 std::vector<int> & tile_union, int & n_near_tiles);
00746
00747
00748
00750 struct EEBriefJet {
00751 double NN_dist;
00752 double kt2;
00753 EEBriefJet * NN;
00754 int _jets_index;
00755
00756 double nx, ny, nz;
00757 };
00758
00760
00761
00762
00764 void _simple_N2_cluster_BriefJet();
00766 void _simple_N2_cluster_EEBriefJet();
00767 };
00768
00769
00770
00771
00772
00773
00774
00775
00776
00777 template<class L> void ClusterSequence::_transfer_input_jets(
00778 const std::vector<L> & pseudojets) {
00779
00780
00781
00782 _jets.reserve(pseudojets.size()*2);
00783
00784
00785
00786
00787
00788 for (unsigned int i = 0; i < pseudojets.size(); i++) {
00789 _jets.push_back(pseudojets[i]);}
00790
00791 }
00792
00793
00794
00795
00796 template<class L> ClusterSequence::ClusterSequence (
00797 const std::vector<L> & pseudojets,
00798 const double & R,
00799 const Strategy & strategy,
00800 const bool & writeout_combinations) {
00801
00802
00803 _transfer_input_jets(pseudojets);
00804
00805
00806 _initialise_and_run(R,strategy,writeout_combinations);
00807 }
00808
00809
00810
00813 template<class L> ClusterSequence::ClusterSequence (
00814 const std::vector<L> & pseudojets,
00815 const JetDefinition & jet_def,
00816 const bool & writeout_combinations) {
00817
00818
00819 _transfer_input_jets(pseudojets);
00820
00821
00822 _initialise_and_run(jet_def,writeout_combinations);
00823 }
00824
00825
00826 inline const std::vector<PseudoJet> & ClusterSequence::jets () const {
00827 return _jets;
00828 }
00829
00830 inline const std::vector<ClusterSequence::history_element> & ClusterSequence::history () const {
00831 return _history;
00832 }
00833
00834 inline unsigned int ClusterSequence::n_particles() const {return _initial_n;}
00835
00836
00837
00838
00839 template <class J> inline void ClusterSequence::_bj_set_jetinfo(
00840 J * const jetA, const int _jets_index) const {
00841 jetA->eta = _jets[_jets_index].rap();
00842 jetA->phi = _jets[_jets_index].phi_02pi();
00843 jetA->kt2 = jet_scale_for_algorithm(_jets[_jets_index]);
00844 jetA->_jets_index = _jets_index;
00845
00846 jetA->NN_dist = _R2;
00847 jetA->NN = NULL;
00848 }
00849
00850
00851
00852
00853
00854 template <class J> inline double ClusterSequence::_bj_dist(
00855 const J * const jetA, const J * const jetB) const {
00856 double dphi = std::abs(jetA->phi - jetB->phi);
00857 double deta = (jetA->eta - jetB->eta);
00858 if (dphi > pi) {dphi = twopi - dphi;}
00859 return dphi*dphi + deta*deta;
00860 }
00861
00862
00863 template <class J> inline double ClusterSequence::_bj_diJ(const J * const jet) const {
00864 double kt2 = jet->kt2;
00865 if (jet->NN != NULL) {if (jet->NN->kt2 < kt2) {kt2 = jet->NN->kt2;}}
00866 return jet->NN_dist * kt2;
00867 }
00868
00869
00870
00871
00872
00873 template <class J> inline void ClusterSequence::_bj_set_NN_nocross(
00874 J * const jet, J * const head, const J * const tail) const {
00875 double NN_dist = _R2;
00876 J * NN = NULL;
00877 if (head < jet) {
00878 for (J * jetB = head; jetB != jet; jetB++) {
00879 double dist = _bj_dist(jet,jetB);
00880 if (dist < NN_dist) {
00881 NN_dist = dist;
00882 NN = jetB;
00883 }
00884 }
00885 }
00886 if (tail > jet) {
00887 for (J * jetB = jet+1; jetB != tail; jetB++) {
00888 double dist = _bj_dist(jet,jetB);
00889 if (dist < NN_dist) {
00890 NN_dist = dist;
00891 NN = jetB;
00892 }
00893 }
00894 }
00895 jet->NN = NN;
00896 jet->NN_dist = NN_dist;
00897 }
00898
00899
00900
00901 template <class J> inline void ClusterSequence::_bj_set_NN_crosscheck(J * const jet,
00902 J * const head, const J * const tail) const {
00903 double NN_dist = _R2;
00904 J * NN = NULL;
00905 for (J * jetB = head; jetB != tail; jetB++) {
00906 double dist = _bj_dist(jet,jetB);
00907 if (dist < NN_dist) {
00908 NN_dist = dist;
00909 NN = jetB;
00910 }
00911 if (dist < jetB->NN_dist) {
00912 jetB->NN_dist = dist;
00913 jetB->NN = jet;
00914 }
00915 }
00916 jet->NN = NN;
00917 jet->NN_dist = NN_dist;
00918 }
00919
00920
00921
00922 FASTJET_END_NAMESPACE
00923
00924 #endif // __FASTJET_CLUSTERSEQUENCE_HH__