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<cmath>
00058 #include "fastjet/Error.hh"
00059 #include "fastjet/JetDefinition.hh"
00060
00061 FASTJET_BEGIN_NAMESPACE
00062
00063
00065 class ClusterSequence {
00066
00067
00068 public:
00069
00071 ClusterSequence () {}
00072
00082 template<class L> ClusterSequence (const std::vector<L> & pseudojets,
00083 const double & R = 1.0,
00084 const Strategy & strategy = Best,
00085 const bool & writeout_combinations = false);
00086
00087
00091 template<class L> ClusterSequence (
00092 const std::vector<L> & pseudojets,
00093 const JetDefinition & jet_def,
00094 const bool & writeout_combinations = false);
00095
00096
00097
00098
00099
00100
00104 std::vector<PseudoJet> inclusive_jets (const double & ptmin = 0.0) const;
00105
00109 int n_exclusive_jets (const double & dcut) const;
00110
00114 std::vector<PseudoJet> exclusive_jets (const double & dcut) const;
00115
00118 std::vector<PseudoJet> exclusive_jets (const int & njets) const;
00119
00122 double exclusive_dmerge (const int & njets) const;
00123
00128 double exclusive_dmerge_max (const int & njets) const;
00129
00140 bool object_in_jet(const PseudoJet & object, const PseudoJet & jet) const;
00141
00147 bool has_parents(const PseudoJet & jet, PseudoJet & parent1,
00148 PseudoJet & parent2) const;
00149
00152 bool has_child(const PseudoJet & jet, PseudoJet & child) const;
00153
00156 bool has_child(const PseudoJet & jet, const PseudoJet * & childp) const;
00157
00161 bool has_partner(const PseudoJet & jet, PseudoJet & partner) const;
00162
00164 std::vector<PseudoJet> constituents (const PseudoJet & jet) const;
00165
00166
00175 void print_jets_for_root(const std::vector<PseudoJet> & jets,
00176 std::ostream & ostr = std::cout) const;
00177
00178
00179
00180
00181
00183 void add_constituents (const PseudoJet & jet,
00184 std::vector<PseudoJet> & subjet_vector) const;
00185
00187 inline Strategy strategy_used () const {return _strategy;}
00188 std::string strategy_string () const;
00189
00191 const JetDefinition & jet_def() const {return _jet_def;}
00192
00196 double jet_scale_for_algorithm(const PseudoJet & jet) const;
00197
00198
00199
00200
00206 void plugin_record_ij_recombination(int jet_i, int jet_j, double dij,
00207 int & newjet_k) {
00208 assert(plugin_activated());
00209 _do_ij_recombination_step(jet_i, jet_j, dij, newjet_k);
00210 }
00211
00215 void plugin_record_ij_recombination(int jet_i, int jet_j, double dij,
00216 const PseudoJet & newjet,
00217 int & newjet_k);
00218
00223 void plugin_record_iB_recombination(int jet_i, double diB) {
00224 assert(plugin_activated());
00225 _do_iB_recombination_step(jet_i, diB);
00226 }
00227
00231 class Extras {
00232 public:
00233 virtual ~Extras() {}
00234 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";}
00235 };
00236
00239 inline void plugin_associate_extras(std::auto_ptr<Extras> extras_in) {
00240 _extras = extras_in;
00241 }
00242
00244 inline bool plugin_activated() const {return _plugin_activated;}
00245
00247 const Extras * extras() const {return _extras.get();}
00248
00249 public:
00253 static void set_jet_algorithm (JetAlgorithm jet_algorithm) {_default_jet_algorithm = jet_algorithm;}
00255 static void set_jet_finder (JetAlgorithm jet_algorithm) {_default_jet_algorithm = jet_algorithm;}
00256
00257
00260 struct history_element{
00261 int parent1;
00262
00263
00264
00265 int parent2;
00266
00267
00268
00269
00270
00271 int child;
00272
00273
00274
00275
00276 int jetp_index;
00277
00278
00279
00280
00281
00282
00283 double dij;
00284
00285
00286 double max_dij_so_far;
00287
00288 };
00289
00290 enum JetType {Invalid=-3, InexistentParent = -2, BeamJet = -1};
00291
00297 const std::vector<PseudoJet> & jets() const;
00298
00301 const std::vector<history_element> & history() const;
00302
00307 unsigned int n_particles() const;
00308
00313 std::vector<int> particle_jet_indices(const std::vector<PseudoJet> &) const;
00314
00330 std::vector<int> unique_history_order() const;
00331
00335 std::vector<PseudoJet> unclustered_particles() const;
00336
00340 void transfer_from_sequence(ClusterSequence & from_seq);
00341
00342 protected:
00343 static JetAlgorithm _default_jet_algorithm;
00344 JetDefinition _jet_def;
00345
00348 bool _potentially_valid(const PseudoJet & jet) const {
00349 return jet.cluster_hist_index() >= 0
00350 && jet.cluster_hist_index() < int(_history.size());
00351 }
00352
00355 template<class L> void _transfer_input_jets(
00356 const std::vector<L> & pseudojets);
00357
00361 void _initialise_and_run (const JetDefinition & jet_def,
00362 const bool & writeout_combinations);
00363
00367 void _initialise_and_run (const double & R,
00368 const Strategy & strategy,
00369 const bool & writeout_combinations);
00370
00373 void _decant_options(const JetDefinition & jet_def,
00374 const bool & writeout_combinations);
00375
00379 void _fill_initial_history();
00380
00384 void _do_ij_recombination_step(const int & jet_i, const int & jet_j,
00385 const double & dij, int & newjet_k);
00386
00389 void _do_iB_recombination_step(const int & jet_i, const double & diB);
00390
00391
00395 std::vector<PseudoJet> _jets;
00396
00397
00401 std::vector<history_element> _history;
00402
00403 bool _writeout_combinations;
00404 int _initial_n;
00405 double _Rparam, _R2, _invR2;
00406 Strategy _strategy;
00407 JetAlgorithm _jet_algorithm;
00408
00409 private:
00410
00411 bool _plugin_activated;
00412 std::auto_ptr<Extras> _extras;
00413
00414 void _really_dumb_cluster ();
00415 void _delaunay_cluster ();
00416 void _simple_N2_cluster ();
00417 void _tiled_N2_cluster ();
00418 void _faster_tiled_N2_cluster ();
00419
00420
00421 void _minheap_faster_tiled_N2_cluster();
00422
00423
00424
00425 void _CP2DChan_cluster();
00426 void _CP2DChan_cluster_2pi2R ();
00427 void _CP2DChan_cluster_2piMultD ();
00428 void _CP2DChan_limited_cluster(double D);
00429 void _do_Cambridge_inclusive_jets();
00430
00431 void _add_step_to_history(const int & step_number, const int & parent1,
00432 const int & parent2, const int & jetp_index,
00433 const double & dij);
00434
00437 void _extract_tree_children(int pos, std::valarray<bool> &,
00438 const std::valarray<int> &, std::vector<int> &) const;
00439
00442 void _extract_tree_parents (int pos, std::valarray<bool> &,
00443 const std::valarray<int> &, std::vector<int> &) const;
00444
00445
00446
00447 typedef std::pair<int,int> TwoVertices;
00448 typedef std::pair<double,TwoVertices> DijEntry;
00449 typedef std::multimap<double,TwoVertices> DistMap;
00450
00452 void _add_ktdistance_to_map(const int & ii,
00453 DistMap & DijMap,
00454 const DynamicNearestNeighbours * DNN);
00455
00457 void _print_banner();
00459 static bool _first_time;
00460
00464 static int _n_exclusive_warnings;
00465
00466
00471 struct BriefJet {
00472 double eta, phi, kt2, NN_dist;
00473 BriefJet * NN;
00474 int _jets_index;
00475 };
00478 class TiledJet {
00479 public:
00480 double eta, phi, kt2, NN_dist;
00481 TiledJet * NN, *previous, * next;
00482 int _jets_index, tile_index, diJ_posn;
00483
00484
00485
00486
00487 inline void label_minheap_update_needed() {diJ_posn = 1;}
00488 inline void label_minheap_update_done() {diJ_posn = 0;}
00489 inline bool minheap_update_needed() const {return diJ_posn==1;}
00490 };
00491
00492
00493
00494
00497 template <class J> void _bj_set_jetinfo( J * const jet,
00498 const int _jets_index) const;
00499
00502 void _bj_remove_from_tiles( TiledJet * const jet) const;
00503
00505 template <class J> double _bj_dist(const J * const jeta,
00506 const J * const jetb) const;
00507
00508
00509
00510 template <class J> double _bj_diJ(const J * const jeta) const;
00511
00515 template <class J> inline J * _bj_of_hindex(
00516 const int hist_index,
00517 J * const head, J * const tail)
00518 const {
00519 J * res;
00520 for(res = head; res<tail; res++) {
00521 if (_jets[res->_jets_index].cluster_hist_index() == hist_index) {break;}
00522 }
00523 return res;
00524 }
00525
00526
00527
00528
00529
00534 template <class J> void _bj_set_NN_nocross(J * const jeta,
00535 J * const head, const J * const tail) const;
00536
00541 template <class J> void _bj_set_NN_crosscheck(J * const jeta,
00542 J * const head, const J * const tail) const;
00543
00544
00545
00548 static const int n_tile_neighbours = 9;
00549
00552 struct Tile {
00554 Tile * begin_tiles[n_tile_neighbours];
00556 Tile ** surrounding_tiles;
00558 Tile ** RH_tiles;
00560 Tile ** end_tiles;
00562 TiledJet * head;
00564 bool tagged;
00565 };
00566 std::vector<Tile> _tiles;
00567 double _tiles_eta_min, _tiles_eta_max;
00568 double _tile_size_eta, _tile_size_phi;
00569 int _n_tiles_phi,_tiles_ieta_min,_tiles_ieta_max;
00570
00571
00572
00573 inline int _tile_index (int ieta, int iphi) const {
00574
00575
00576 return (ieta-_tiles_ieta_min)*_n_tiles_phi
00577 + (iphi+_n_tiles_phi) % _n_tiles_phi;
00578 }
00579
00580
00581
00582 int _tile_index(const double & eta, const double & phi) const;
00583 void _tj_set_jetinfo ( TiledJet * const jet, const int _jets_index);
00584 void _bj_remove_from_tiles(TiledJet * const jet);
00585 void _initialise_tiles();
00586 void _print_tiles(TiledJet * briefjets ) const;
00587 void _add_neighbours_to_tile_union(const int tile_index,
00588 std::vector<int> & tile_union, int & n_near_tiles) const;
00589 void _add_untagged_neighbours_to_tile_union(const int tile_index,
00590 std::vector<int> & tile_union, int & n_near_tiles);
00591
00592
00593 };
00594
00595
00596
00597
00598
00599
00600
00601
00602
00603
00604 template<class L> void ClusterSequence::_transfer_input_jets(
00605 const std::vector<L> & pseudojets) {
00606
00607
00608
00609 _jets.reserve(pseudojets.size()*2);
00610
00611
00612
00613
00614
00615 for (unsigned int i = 0; i < pseudojets.size(); i++) {
00616 _jets.push_back(pseudojets[i]);}
00617
00618 }
00619
00620
00621
00622
00623 template<class L> ClusterSequence::ClusterSequence (
00624 const std::vector<L> & pseudojets,
00625 const double & R,
00626 const Strategy & strategy,
00627 const bool & writeout_combinations) {
00628
00629
00630 _transfer_input_jets(pseudojets);
00631
00632
00633 _initialise_and_run(R,strategy,writeout_combinations);
00634 }
00635
00636
00637
00640 template<class L> ClusterSequence::ClusterSequence (
00641 const std::vector<L> & pseudojets,
00642 const JetDefinition & jet_def,
00643 const bool & writeout_combinations) {
00644
00645
00646 _transfer_input_jets(pseudojets);
00647
00648
00649 _initialise_and_run(jet_def,writeout_combinations);
00650 }
00651
00652
00653 inline const std::vector<PseudoJet> & ClusterSequence::jets () const {
00654 return _jets;
00655 }
00656
00657 inline const std::vector<ClusterSequence::history_element> & ClusterSequence::history () const {
00658 return _history;
00659 }
00660
00661 inline unsigned int ClusterSequence::n_particles() const {return _initial_n;}
00662
00663
00664
00665
00666 template <class J> inline void ClusterSequence::_bj_set_jetinfo(
00667 J * const jetA, const int _jets_index) const {
00668 jetA->eta = _jets[_jets_index].rap();
00669 jetA->phi = _jets[_jets_index].phi_02pi();
00670 jetA->kt2 = jet_scale_for_algorithm(_jets[_jets_index]);
00671 jetA->_jets_index = _jets_index;
00672
00673 jetA->NN_dist = _R2;
00674 jetA->NN = NULL;
00675 }
00676
00677
00678
00679
00680
00681 template <class J> inline double ClusterSequence::_bj_dist(
00682 const J * const jetA, const J * const jetB) const {
00683 double dphi = std::abs(jetA->phi - jetB->phi);
00684 double deta = (jetA->eta - jetB->eta);
00685 if (dphi > pi) {dphi = twopi - dphi;}
00686 return dphi*dphi + deta*deta;
00687 }
00688
00689
00690 template <class J> inline double ClusterSequence::_bj_diJ(const J * const jet) const {
00691 double kt2 = jet->kt2;
00692 if (jet->NN != NULL) {if (jet->NN->kt2 < kt2) {kt2 = jet->NN->kt2;}}
00693 return jet->NN_dist * kt2;
00694 }
00695
00696
00697
00698
00699
00700 template <class J> inline void ClusterSequence::_bj_set_NN_nocross(
00701 J * const jet, J * const head, const J * const tail) const {
00702 double NN_dist = _R2;
00703 J * NN = NULL;
00704 if (head < jet) {
00705 for (J * jetB = head; jetB != jet; jetB++) {
00706 double dist = _bj_dist(jet,jetB);
00707 if (dist < NN_dist) {
00708 NN_dist = dist;
00709 NN = jetB;
00710 }
00711 }
00712 }
00713 if (tail > jet) {
00714 for (J * jetB = jet+1; jetB != tail; jetB++) {
00715 double dist = _bj_dist(jet,jetB);
00716 if (dist < NN_dist) {
00717 NN_dist = dist;
00718 NN = jetB;
00719 }
00720 }
00721 }
00722 jet->NN = NN;
00723 jet->NN_dist = NN_dist;
00724 }
00725
00726
00727
00728 template <class J> inline void ClusterSequence::_bj_set_NN_crosscheck(J * const jet,
00729 J * const head, const J * const tail) const {
00730 double NN_dist = _R2;
00731 J * NN = NULL;
00732 for (J * jetB = head; jetB != tail; jetB++) {
00733 double dist = _bj_dist(jet,jetB);
00734 if (dist < NN_dist) {
00735 NN_dist = dist;
00736 NN = jetB;
00737 }
00738 if (dist < jetB->NN_dist) {
00739 jetB->NN_dist = dist;
00740 jetB->NN = jet;
00741 }
00742 }
00743 jet->NN = NN;
00744 jet->NN_dist = NN_dist;
00745 }
00746
00747
00748
00749
00750 FASTJET_END_NAMESPACE
00751
00752 #endif // __FASTJET_CLUSTERSEQUENCE_HH__