fastjet/internal/LazyTiling9.hh

0001 #ifndef __FASTJET_LAZYTILING9_HH__
0002 #define __FASTJET_LAZYTILING9_HH__
0003
0004 //#define INSTRUMENT2 1
0005
0006 //FJSTARTHEADER
0007 // $Id$
0008 //
0009 // Copyright (c) 2005-2021, Matteo Cacciari, Gavin P. Salam and Gregory Soyez
0010 //
0011 //----------------------------------------------------------------------
0012 // This file is part of FastJet.
0013 //
0014 //  FastJet is free software; you can redistribute it and/or modify
0015 //  it under the terms of the GNU General Public License as published by
0016 //  the Free Software Foundation; either version 2 of the License, or
0017 //  (at your option) any later version.
0018 //
0019 //  The algorithms that underlie FastJet have required considerable
0020 //  development. They are described in the original FastJet paper,
0021 //  hep-ph/0512210 and in the manual, arXiv:1111.6097. If you use
0022 //  FastJet as part of work towards a scientific publication, please
0023 //  quote the version you use and include a citation to the manual and
0024 //  optionally also to hep-ph/0512210.
0025 //
0026 //  FastJet is distributed in the hope that it will be useful,
0027 //  but WITHOUT ANY WARRANTY; without even the implied warranty of
0028 //  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
0029 //  GNU General Public License for more details.
0030 //
0031 //  You should have received a copy of the GNU General Public License
0032 //  along with FastJet. If not, see <http://www.gnu.org/licenses/>.
0033 //----------------------------------------------------------------------
0034 //FJENDHEADER
0035
0036 //#include "fastjet/PseudoJet.hh"
0037 #include "fastjet/internal/MinHeap.hh"
0038 #include "fastjet/ClusterSequence.hh"
0039 #include "fastjet/internal/LazyTiling9Alt.hh"
0040
0041
0042 FASTJET_BEGIN_NAMESPACE      // defined in fastjet/internal/base.hh
0043 template<int NN>
0044 class Tile2Base {
0045 public:
0046   /// pointers to neighbouring tiles, including self
0047   Tile2Base *   begin_tiles[NN];
0048   /// neighbouring tiles, excluding self
0049   Tile2Base **  surrounding_tiles;
0050   /// half of neighbouring tiles, no self
0051   Tile2Base **  RH_tiles;
0052   /// just beyond end of tiles
0053   Tile2Base **  end_tiles;
0054   /// start of list of BriefJets contained in this tile
0055   TiledJet * head;
0056   /// sometimes useful to be able to tag a tile
0057   bool     tagged;
0058   /// true for tiles where the delta phi calculation needs
0059   /// potentially to account for periodicity in phi
0060   bool     use_periodic_delta_phi;
0061   /// for all particles in the tile, this stores the largest of the
0062   /// (squared) nearest-neighbour distances.
0063   double max_NN_dist;
0064   double eta_centre, phi_centre;
0065
0066   /// returns the number of jets in the tile; useful principally for
0067   /// diagnostics
0068   int jet_count() const {
0069     int count = 0;
0070     const TiledJet * jet = head;
0071     while (jet != 0) {
0072       count++;
0073       jet = jet->next;
0074     }
0075     return count;
0076   }
0077 };
0078
0079
0080 typedef Tile2Base<9> Tile2;
0081
0082 // class Tile2 : public Tile2Base {
0083 // public:
0084 //   /// pointers to neighbouring tiles, including self
0085 //   Tile2 *   begin_tiles[n_tile_neighbours];
0086 //   bool is_near_zero_phi(double tile_size_phi) const {
0087 //     return phi_centre < tile_size_phi || (twopi-phi_centre) < tile_size_phi;
0088 //   }
0089 // };
0090
0091 // class Tile2  {
0092 // public:
0093 //   /// pointers to neighbouring tiles, including self
0094 //   Tile2 *   begin_tiles[n_tile_neighbours];
0095 //   /// neighbouring tiles, excluding self
0096 //   Tile2 **  surrounding_tiles;
0097 //   /// half of neighbouring tiles, no self
0098 //   Tile2 **  RH_tiles;
0099 //   /// just beyond end of tiles
0100 //   Tile2 **  end_tiles;
0101 //   /// start of list of BriefJets contained in this tile
0102 //   TiledJet * head;
0103 //   /// sometimes useful to be able to tag a tile
0104 //   bool     tagged;
0105 //   /// for all particles in the tile, this stores the largest of the
0106 //   /// (squared) nearest-neighbour distances.
0107 //   double max_NN_dist;
0108 //   double eta_centre, phi_centre;
0109 //   bool is_near_zero_phi(double tile_size_phi) const {
0110 //     return phi_centre < tile_size_phi || (twopi-phi_centre) < tile_size_phi;
0111 //   }
0112 // };
0113
0114
0115
0116 //----------------------------------------------------------------------
0117 class LazyTiling9 {
0118 public:
0119   LazyTiling9(ClusterSequence & cs);
0120
0121   void run();
0122
0123   //void get_next_clustering(int & jetA_index, int & jetB_index, double & dij);
0124
0125
0126 protected:
0127   ClusterSequence & _cs;
0128   const std::vector<PseudoJet> & _jets;
0129   std::vector<Tile2> _tiles;
0130
0131 #ifdef INSTRUMENT2
0132   int _ncall; // GPS tmp
0133   int _ncall_dtt; // GPS tmp
0134 #endif // INSTRUMENT2
0135
0136   double _Rparam, _R2, _invR2;
0137   double _tiles_eta_min, _tiles_eta_max;
0138   double _tile_size_eta, _tile_size_phi;
0139   double _tile_half_size_eta, _tile_half_size_phi;
0140   int    _n_tiles_phi,_tiles_ieta_min,_tiles_ieta_max;
0141
0142   std::vector<TiledJet *> _jets_for_minheap;
0143
0144   //MinHeap _minheap;
0145
0146   void _initialise_tiles();
0147
0148   // reasonably robust return of tile index given ieta and iphi, in particular
0149   // it works even if iphi is negative
0150   inline int _tile_index (int ieta, int iphi) const {
0151     // note that (-1)%n = -1 so that we have to add _n_tiles_phi
0152     // before performing modulo operation
0153     return (ieta-_tiles_ieta_min)*_n_tiles_phi
0154                   + (iphi+_n_tiles_phi) % _n_tiles_phi;
0155   }
0156
0157   void  _bj_remove_from_tiles(TiledJet * const jet);
0158
0159   /// returns the tile index given the eta and phi values of a jet
0160   int _tile_index(const double eta, const double phi) const;
0161
0162   // sets up information regarding the tiling of the given jet
0163   void _tj_set_jetinfo(TiledJet * const jet, const int _jets_index);
0164
0165   void _print_tiles(TiledJet * briefjets ) const;
0166   void _add_neighbours_to_tile_union(const int tile_index,
0167          std::vector<int> & tile_union, int & n_near_tiles) const;
0168   void _add_untagged_neighbours_to_tile_union(const int tile_index,
0169          std::vector<int> & tile_union, int & n_near_tiles);
0170   void _add_untagged_neighbours_to_tile_union_using_max_info(const TiledJet * const jet,
0171          std::vector<int> & tile_union, int & n_near_tiles);
0172   double _distance_to_tile(const TiledJet * bj, const Tile2 *)
0173 #ifdef INSTRUMENT2
0174     ;
0175 #else
0176     const;
0177 #endif
0178   void _update_jetX_jetI_NN(TiledJet * jetX, TiledJet * jetI, std::vector<TiledJet *> & jets_for_minheap);
0179
0180   void _set_NN(TiledJet * jetI, std::vector<TiledJet *> & jets_for_minheap);
0181
0182   // return the diJ (multiplied by _R2) for this jet assuming its NN
0183   // info is correct
0184   template <class J> double _bj_diJ(const J * const jet) const {
0185     double kt2 = jet->kt2;
0186     if (jet->NN != NULL) {if (jet->NN->kt2 < kt2) {kt2 = jet->NN->kt2;}}
0187     return jet->NN_dist * kt2;
0188   }
0189
0190
0191   //----------------------------------------------------------------------
0192   template <class J> inline void _bj_set_jetinfo(
0193                             J * const jetA, const int _jets_index) const {
0194     jetA->eta  = _jets[_jets_index].rap();
0195     jetA->phi  = _jets[_jets_index].phi_02pi();
0196     jetA->kt2  = _cs.jet_scale_for_algorithm(_jets[_jets_index]);
0197     jetA->_jets_index = _jets_index;
0198     // initialise NN info as well
0199     jetA->NN_dist = _R2;
0200     jetA->NN      = NULL;
0201   }
0202
0203
0204   //----------------------------------------------------------------------
0205   template <class J> inline double _bj_dist(
0206                 const J * const jetA, const J * const jetB)
0207 #ifdef INSTRUMENT2
0208     {
0209     _ncall++; // GPS tmp
0210 #else
0211     const {
0212 #endif
0213     double dphi = std::abs(jetA->phi - jetB->phi);
0214     double deta = (jetA->eta - jetB->eta);
0215     if (dphi > pi) {dphi = twopi - dphi;}
0216     return dphi*dphi + deta*deta;
0217   }
0218
0219
0220   //----------------------------------------------------------------------
0221   template <class J> inline double _bj_dist_not_periodic(
0222                 const J * const jetA, const J * const jetB)
0223 #ifdef INSTRUMENT2
0224     {
0225     _ncall++; // GPS tmp
0226 #else
0227     const {
0228 #endif
0229     //_ncall++; // GPS tmp
0230     double dphi = jetA->phi - jetB->phi;
0231     double deta = (jetA->eta - jetB->eta);
0232     return dphi*dphi + deta*deta;
0233   }
0234
0235 };
0236
0237
0238 FASTJET_END_NAMESPACE
0239
0240 #endif // __FASTJET_LAZYTILING9_HH__