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