1 #ifndef __FASTJET_LAZYTILING9SEPARATEGHOSTS_HH__ 2 #define __FASTJET_LAZYTILING9SEPARATEGHOSTS_HH__ 3 4 //FJSTARTHEADER 5 // $Id: LazyTiling9SeparateGhosts.hh 4442 2020-05-05 07:50:11Z soyez $ 6 // 7 // Copyright (c) 2005-2020, Matteo Cacciari, Gavin P. Salam and Gregory Soyez 8 // 9 //---------------------------------------------------------------------- 10 // This file is part of FastJet. 11 // 12 // FastJet is free software; you can redistribute it and/or modify 13 // it under the terms of the GNU General Public License as published by 14 // the Free Software Foundation; either version 2 of the License, or 15 // (at your option) any later version. 16 // 17 // The algorithms that underlie FastJet have required considerable 18 // development. They are described in the original FastJet paper, 19 // hep-ph/0512210 and in the manual, arXiv:1111.6097. If you use 20 // FastJet as part of work towards a scientific publication, please 21 // quote the version you use and include a citation to the manual and 22 // optionally also to hep-ph/0512210. 23 // 24 // FastJet is distributed in the hope that it will be useful, 25 // but WITHOUT ANY WARRANTY; without even the implied warranty of 26 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the 27 // GNU General Public License for more details. 28 // 29 // You should have received a copy of the GNU General Public License 30 // along with FastJet. If not, see <http://www.gnu.org/licenses/>. 31 //---------------------------------------------------------------------- 32 //FJENDHEADER 33 34 //#include "fastjet/PseudoJet.hh" 35 #include "fastjet/internal/MinHeap.hh" 36 #include "fastjet/ClusterSequence.hh" 37 #include "fastjet/internal/LazyTiling9Alt.hh" 38 39 FASTJET_BEGIN_NAMESPACE // defined in fastjet/internal/base.hh 40 41 class TiledJet3 { 42 public: 43 double eta, phi, kt2, NN_dist; 44 TiledJet3 * NN, *previous, * next; 45 int _jets_index, tile_index; 46 bool _minheap_update_needed; 47 bool is_ghost; 48 49 // indicate whether jets need to have their minheap entries 50 // updated). label_minheap_update_needed()51 inline void label_minheap_update_needed() {_minheap_update_needed = true;} label_minheap_update_done()52 inline void label_minheap_update_done() {_minheap_update_needed = false;} minheap_update_needed() const53 inline bool minheap_update_needed() const {return _minheap_update_needed;} 54 }; 55 56 57 class Tile3 { 58 public: 59 /// pointers to neighbouring tiles, including self 60 Tile3 * begin_tiles[n_tile_neighbours]; 61 /// neighbouring tiles, excluding self 62 Tile3 ** surrounding_tiles; 63 /// half of neighbouring tiles, no self 64 Tile3 ** RH_tiles; 65 /// just beyond end of tiles 66 Tile3 ** end_tiles; 67 /// start of list of BriefJets contained in this tile 68 TiledJet3 * head; 69 /// start of list of BriefJets contained in this tile 70 TiledJet3 * ghost_head; 71 /// sometimes useful to be able to tag a tile 72 bool tagged; 73 /// for all particles in the tile, this stores the largest of the 74 /// (squared) nearest-neighbour distances. 75 double max_NN_dist; 76 double eta_centre, phi_centre; 77 is_near_zero_phi(double tile_size_phi) const78 bool is_near_zero_phi(double tile_size_phi) const { 79 return phi_centre < tile_size_phi || (twopi-phi_centre) < tile_size_phi; 80 } 81 }; 82 83 84 //---------------------------------------------------------------------- 85 class LazyTiling9SeparateGhosts { 86 public: 87 LazyTiling9SeparateGhosts(ClusterSequence & cs); 88 89 void run(); 90 91 //void get_next_clustering(int & jetA_index, int & jetB_index, double & dij); 92 93 /// this is the pt2 threshold below which particles will be 94 /// considered as ghosts. 95 /// 96 /// Note that as it stands, a user can decide to change that 97 /// value. Note however that this has to be done at the user's own 98 /// risk (this is an internal part of fastjet).In a similar spirit, 99 /// the interface to access this valus might also change in a future 100 /// release of FastJet. 101 static double ghost_pt2_threshold; 102 103 protected: 104 ClusterSequence & _cs; 105 const std::vector<PseudoJet> & _jets; 106 std::vector<Tile3> _tiles; 107 108 109 double _Rparam, _R2, _invR2; 110 double _tiles_eta_min, _tiles_eta_max; 111 double _tile_size_eta, _tile_size_phi; 112 double _tile_half_size_eta, _tile_half_size_phi; 113 int _n_tiles_phi,_tiles_ieta_min,_tiles_ieta_max; 114 115 std::vector<TiledJet3 *> _jets_for_minheap; 116 117 //MinHeap _minheap; 118 119 void _initialise_tiles(); 120 121 // reasonably robust return of tile index given ieta and iphi, in particular 122 // it works even if iphi is negative _tile_index(int ieta,int iphi) const123 inline int _tile_index (int ieta, int iphi) const { 124 // note that (-1)%n = -1 so that we have to add _n_tiles_phi 125 // before performing modulo operation 126 return (ieta-_tiles_ieta_min)*_n_tiles_phi 127 + (iphi+_n_tiles_phi) % _n_tiles_phi; 128 } 129 130 void _bj_remove_from_tiles(TiledJet3 * const jet); 131 132 /// returns the tile index given the eta and phi values of a jet 133 int _tile_index(const double eta, const double phi) const; 134 135 // sets up information regarding the tiling of the given jet 136 void _tj_set_jetinfo(TiledJet3 * const jet, const int _jets_index, bool is_ghost); 137 138 void _print_tiles(TiledJet3 * briefjets ) const; 139 void _add_neighbours_to_tile_union(const int tile_index, 140 std::vector<int> & tile_union, int & n_near_tiles) const; 141 void _add_untagged_neighbours_to_tile_union(const int tile_index, 142 std::vector<int> & tile_union, int & n_near_tiles); 143 void _add_untagged_neighbours_to_tile_union_using_max_info(const TiledJet3 * const jet, 144 std::vector<int> & tile_union, int & n_near_tiles); 145 double _distance_to_tile(const TiledJet3 * bj, const Tile3 *) const; 146 void _update_jetX_jetI_NN(TiledJet3 * jetX, TiledJet3 * jetI, std::vector<TiledJet3 *> & jets_for_minheap); 147 148 void _set_NN(TiledJet3 * jetI, std::vector<TiledJet3 *> & jets_for_minheap); 149 150 // return the diJ (multiplied by _R2) for this jet assuming its NN 151 // info is correct _bj_diJ(const J * const jet) const152 template <class J> double _bj_diJ(const J * const jet) const { 153 double kt2 = jet->kt2; 154 if (jet->NN != NULL) {if (jet->NN->kt2 < kt2) {kt2 = jet->NN->kt2;}} 155 return jet->NN_dist * kt2; 156 } 157 158 159 //---------------------------------------------------------------------- _bj_set_jetinfo(J * const jetA,const int _jets_index) const160 template <class J> inline void _bj_set_jetinfo( 161 J * const jetA, const int _jets_index) const { 162 jetA->eta = _jets[_jets_index].rap(); 163 jetA->phi = _jets[_jets_index].phi_02pi(); 164 jetA->kt2 = _cs.jet_scale_for_algorithm(_jets[_jets_index]); 165 jetA->_jets_index = _jets_index; 166 // initialise NN info as well 167 jetA->NN_dist = _R2; 168 jetA->NN = NULL; 169 } 170 171 172 //---------------------------------------------------------------------- _bj_dist(const J * const jetA,const J * const jetB) const173 template <class J> inline double _bj_dist( 174 const J * const jetA, const J * const jetB) const { 175 double dphi = std::abs(jetA->phi - jetB->phi); 176 double deta = (jetA->eta - jetB->eta); 177 if (dphi > pi) {dphi = twopi - dphi;} 178 return dphi*dphi + deta*deta; 179 } 180 181 182 //---------------------------------------------------------------------- _bj_dist_not_periodic(const J * const jetA,const J * const jetB) const183 template <class J> inline double _bj_dist_not_periodic( 184 const J * const jetA, const J * const jetB) const { 185 double dphi = jetA->phi - jetB->phi; 186 double deta = (jetA->eta - jetB->eta); 187 return dphi*dphi + deta*deta; 188 } 189 190 }; 191 192 193 FASTJET_END_NAMESPACE 194 195 #endif // __FASTJET_LAZYTILING9SEPARATEGHOSTS_HH__ 196