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
00131 std::vector<PseudoJet> constituents (const PseudoJet & jet) const;
00133 void add_constituents (const PseudoJet & jet,
00134 std::vector<PseudoJet> & subjet_vector) const;
00135
00137 inline Strategy strategy_used () const {return _strategy;}
00138 std::string strategy_string () const;
00139
00140
00144 double jet_scale_for_algorithm(const PseudoJet & jet) const;
00145
00146
00147
00148
00154 void plugin_record_ij_recombination(int jet_i, int jet_j, double dij,
00155 int & newjet_k) {
00156 assert(plugin_activated());
00157 _do_ij_recombination_step(jet_i, jet_j, dij, newjet_k);
00158 }
00159
00163 void plugin_record_ij_recombination(int jet_i, int jet_j, double dij,
00164 const PseudoJet & newjet,
00165 int & newjet_k);
00166
00171 void plugin_record_iB_recombination(int jet_i, double diB) {
00172 assert(plugin_activated());
00173 _do_iB_recombination_step(jet_i, diB);
00174 }
00175
00179 class Extras {
00180 public:
00181 virtual ~Extras() {}
00182 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";}
00183 };
00184
00187 inline void plugin_associate_extras(std::auto_ptr<Extras> extras_in) {
00188 _extras = extras_in;
00189 }
00190
00192 inline bool plugin_activated() const {return _plugin_activated;}
00193
00195 const Extras * extras() const {return _extras.get();}
00196
00197 public:
00201 static void set_jet_finder (JetFinder jet_finder) {_default_jet_finder = jet_finder;}
00202
00203
00206 struct history_element{
00207 int parent1;
00208
00209
00210
00211 int parent2;
00212
00213
00214
00215
00216
00217 int child;
00218
00219
00220
00221
00222 int jetp_index;
00223
00224
00225
00226
00227
00228
00229 double dij;
00230
00231
00232 double max_dij_so_far;
00233
00234 };
00235
00236 enum JetType {Invalid=-3, InexistentParent = -2, BeamJet = -1};
00237
00243 const std::vector<PseudoJet> & jets() const;
00244
00247 const std::vector<history_element> & history() const;
00248
00253 unsigned int n_particles() const;
00254
00259 std::vector<int> particle_jet_indices(const std::vector<PseudoJet> &) const;
00260
00276 std::vector<int> unique_history_order() const;
00277
00281 std::vector<PseudoJet> unclustered_particles() const;
00282
00283
00284 protected:
00285 static JetFinder _default_jet_finder;
00286 JetDefinition _jet_def;
00287
00290 template<class L> void _transfer_input_jets(
00291 const std::vector<L> & pseudojets);
00292
00296 void _initialise_and_run (const JetDefinition & jet_def,
00297 const bool & writeout_combinations);
00298
00302 void _initialise_and_run (const double & R,
00303 const Strategy & strategy,
00304 const bool & writeout_combinations);
00305
00308 void _decant_options(const JetDefinition & jet_def,
00309 const bool & writeout_combinations);
00310
00314 void _fill_initial_history();
00315
00319 void _do_ij_recombination_step(const int & jet_i, const int & jet_j,
00320 const double & dij, int & newjet_k);
00321
00324 void _do_iB_recombination_step(const int & jet_i, const double & diB);
00325
00326
00330 std::vector<PseudoJet> _jets;
00331
00332
00336 std::vector<history_element> _history;
00337
00338 bool _writeout_combinations;
00339 int _initial_n;
00340 double _Rparam, _R2, _invR2;
00341 Strategy _strategy;
00342 JetFinder _jet_finder;
00343
00344 private:
00345
00346 bool _plugin_activated;
00347 std::auto_ptr<Extras> _extras;
00348
00349 void _really_dumb_cluster ();
00350 void _delaunay_cluster ();
00351 void _simple_N2_cluster ();
00352 void _tiled_N2_cluster ();
00353 void _faster_tiled_N2_cluster ();
00354
00355
00356 void _minheap_faster_tiled_N2_cluster();
00357
00358
00359
00360 void _CP2DChan_cluster();
00361 void _CP2DChan_cluster_2pi2R ();
00362 void _CP2DChan_cluster_2piMultD ();
00363 void _CP2DChan_limited_cluster(double D);
00364 void _do_Cambridge_inclusive_jets();
00365
00366 void _add_step_to_history(const int & step_number, const int & parent1,
00367 const int & parent2, const int & jetp_index,
00368 const double & dij);
00369
00372 void _extract_tree_children(int pos, std::valarray<bool> &,
00373 const std::valarray<int> &, std::vector<int> &) const;
00374
00377 void _extract_tree_parents (int pos, std::valarray<bool> &,
00378 const std::valarray<int> &, std::vector<int> &) const;
00379
00380
00381
00382 typedef std::pair<int,int> TwoVertices;
00383 typedef std::pair<double,TwoVertices> DijEntry;
00384 typedef std::multimap<double,TwoVertices> DistMap;
00385
00387 void _add_ktdistance_to_map(const int & ii,
00388 DistMap & DijMap,
00389 const DynamicNearestNeighbours * DNN);
00390
00392 void _print_banner();
00394 static bool _first_time;
00395
00399 static int _n_exclusive_warnings;
00400
00401
00406 struct BriefJet {
00407 double eta, phi, kt2, NN_dist;
00408 BriefJet * NN;
00409 int _jets_index;
00410 };
00413 class TiledJet {
00414 public:
00415 double eta, phi, kt2, NN_dist;
00416 TiledJet * NN, *previous, * next;
00417 int _jets_index, tile_index, diJ_posn;
00418
00419
00420
00421
00422 inline void label_minheap_update_needed() {diJ_posn = 1;}
00423 inline void label_minheap_update_done() {diJ_posn = 0;}
00424 inline bool minheap_update_needed() const {return diJ_posn==1;}
00425 };
00426
00427
00428
00429
00432 template <class J> void _bj_set_jetinfo( J * const jet,
00433 const int _jets_index) const;
00434
00437 void _bj_remove_from_tiles( TiledJet * const jet) const;
00438
00440 template <class J> double _bj_dist(const J * const jeta,
00441 const J * const jetb) const;
00442
00443
00444
00445 template <class J> double _bj_diJ(const J * const jeta) const;
00446
00450 template <class J> inline J * _bj_of_hindex(
00451 const int hist_index,
00452 J * const head, J * const tail)
00453 const {
00454 J * res;
00455 for(res = head; res<tail; res++) {
00456 if (_jets[res->_jets_index].cluster_hist_index() == hist_index) {break;}
00457 }
00458 return res;
00459 }
00460
00461
00462
00463
00464
00469 template <class J> void _bj_set_NN_nocross(J * const jeta,
00470 J * const head, const J * const tail) const;
00471
00476 template <class J> void _bj_set_NN_crosscheck(J * const jeta,
00477 J * const head, const J * const tail) const;
00478
00479
00480
00483 static const int n_tile_neighbours = 9;
00484
00487 struct Tile {
00489 Tile * begin_tiles[n_tile_neighbours];
00491 Tile ** surrounding_tiles;
00493 Tile ** RH_tiles;
00495 Tile ** end_tiles;
00497 TiledJet * head;
00499 bool tagged;
00500 };
00501 std::vector<Tile> _tiles;
00502 double _tiles_eta_min, _tiles_eta_max;
00503 double _tile_size_eta, _tile_size_phi;
00504 int _n_tiles_phi,_tiles_ieta_min,_tiles_ieta_max;
00505
00506
00507
00508 inline int _tile_index (int ieta, int iphi) const {
00509
00510
00511 return (ieta-_tiles_ieta_min)*_n_tiles_phi
00512 + (iphi+_n_tiles_phi) % _n_tiles_phi;
00513 }
00514
00515
00516
00517 int _tile_index(const double & eta, const double & phi) const;
00518 void _tj_set_jetinfo ( TiledJet * const jet, const int _jets_index);
00519 void _bj_remove_from_tiles(TiledJet * const jet);
00520 void _initialise_tiles();
00521 void _print_tiles(TiledJet * briefjets ) const;
00522 void _add_neighbours_to_tile_union(const int tile_index,
00523 std::vector<int> & tile_union, int & n_near_tiles) const;
00524 void _add_untagged_neighbours_to_tile_union(const int tile_index,
00525 std::vector<int> & tile_union, int & n_near_tiles);
00526
00527
00528 };
00529
00530
00531
00532
00533
00534
00535
00536
00537
00538
00539 template<class L> void ClusterSequence::_transfer_input_jets(
00540 const std::vector<L> & pseudojets) {
00541
00542
00543
00544 _jets.reserve(pseudojets.size()*2);
00545
00546
00547
00548
00549
00550 for (unsigned int i = 0; i < pseudojets.size(); i++) {
00551 _jets.push_back(pseudojets[i]);}
00552
00553 }
00554
00555
00556
00557
00558 template<class L> ClusterSequence::ClusterSequence (
00559 const std::vector<L> & pseudojets,
00560 const double & R,
00561 const Strategy & strategy,
00562 const bool & writeout_combinations) {
00563
00564
00565 _transfer_input_jets(pseudojets);
00566
00567
00568 _initialise_and_run(R,strategy,writeout_combinations);
00569 }
00570
00571
00572
00575 template<class L> ClusterSequence::ClusterSequence (
00576 const std::vector<L> & pseudojets,
00577 const JetDefinition & jet_def,
00578 const bool & writeout_combinations) {
00579
00580
00581 _transfer_input_jets(pseudojets);
00582
00583
00584 _initialise_and_run(jet_def,writeout_combinations);
00585 }
00586
00587
00588 inline const std::vector<PseudoJet> & ClusterSequence::jets () const {
00589 return _jets;
00590 }
00591
00592 inline const std::vector<ClusterSequence::history_element> & ClusterSequence::history () const {
00593 return _history;
00594 }
00595
00596 inline unsigned int ClusterSequence::n_particles() const {return _initial_n;}
00597
00598
00599
00600 inline double ClusterSequence::jet_scale_for_algorithm(
00601 const PseudoJet & jet) const {
00602 if (_jet_finder == kt_algorithm) {return jet.kt2();}
00603 else if (_jet_finder == cambridge_algorithm) {return 1.0;}
00604 else {throw Error("Unrecognised jet finder");}
00605 }
00606
00607
00608
00609 template <class J> inline void ClusterSequence::_bj_set_jetinfo(
00610 J * const jetA, const int _jets_index) const {
00611 jetA->eta = _jets[_jets_index].rap();
00612 jetA->phi = _jets[_jets_index].phi_02pi();
00613 jetA->kt2 = jet_scale_for_algorithm(_jets[_jets_index]);
00614 jetA->_jets_index = _jets_index;
00615
00616 jetA->NN_dist = _R2;
00617 jetA->NN = NULL;
00618 }
00619
00620
00621
00622
00623
00624 template <class J> inline double ClusterSequence::_bj_dist(
00625 const J * const jetA, const J * const jetB) const {
00626 double dphi = std::abs(jetA->phi - jetB->phi);
00627 double deta = (jetA->eta - jetB->eta);
00628 if (dphi > pi) {dphi = twopi - dphi;}
00629 return dphi*dphi + deta*deta;
00630 }
00631
00632
00633 template <class J> inline double ClusterSequence::_bj_diJ(const J * const jet) const {
00634 double kt2 = jet->kt2;
00635 if (jet->NN != NULL) {if (jet->NN->kt2 < kt2) {kt2 = jet->NN->kt2;}}
00636 return jet->NN_dist * kt2;
00637 }
00638
00639
00640
00641
00642
00643 template <class J> inline void ClusterSequence::_bj_set_NN_nocross(
00644 J * const jet, J * const head, const J * const tail) const {
00645 double NN_dist = _R2;
00646 J * NN = NULL;
00647 if (head < jet) {
00648 for (J * jetB = head; jetB != jet; jetB++) {
00649 double dist = _bj_dist(jet,jetB);
00650 if (dist < NN_dist) {
00651 NN_dist = dist;
00652 NN = jetB;
00653 }
00654 }
00655 }
00656 if (tail > jet) {
00657 for (J * jetB = jet+1; jetB != tail; jetB++) {
00658 double dist = _bj_dist(jet,jetB);
00659 if (dist < NN_dist) {
00660 NN_dist = dist;
00661 NN = jetB;
00662 }
00663 }
00664 }
00665 jet->NN = NN;
00666 jet->NN_dist = NN_dist;
00667 }
00668
00669
00670
00671 template <class J> inline void ClusterSequence::_bj_set_NN_crosscheck(J * const jet,
00672 J * const head, const J * const tail) const {
00673 double NN_dist = _R2;
00674 J * NN = NULL;
00675 for (J * jetB = head; jetB != tail; jetB++) {
00676 double dist = _bj_dist(jet,jetB);
00677 if (dist < NN_dist) {
00678 NN_dist = dist;
00679 NN = jetB;
00680 }
00681 if (dist < jetB->NN_dist) {
00682 jetB->NN_dist = dist;
00683 jetB->NN = jet;
00684 }
00685 }
00686 jet->NN = NN;
00687 jet->NN_dist = NN_dist;
00688 }
00689
00690
00691
00692
00693 FASTJET_END_NAMESPACE
00694
00695 #endif // __FASTJET_CLUSTERSEQUENCE_HH__