fastjet/internal/LazyTiling9Alt.hh

0001 #ifndef __FASTJET_LAZYTILING9ALT_HH__
0002 #define __FASTJET_LAZYTILING9ALT_HH__
0003
0004 //FJSTARTHEADER
0005 // $Id$
0006 //
0007 // Copyright (c) 2005-2021, Matteo Cacciari, Gavin P. Salam and Gregory Soyez
0008 //
0009 //----------------------------------------------------------------------
0010 // This file is part of FastJet.
0011 //
0012 //  FastJet is free software; you can redistribute it and/or modify
0013 //  it under the terms of the GNU General Public License as published by
0014 //  the Free Software Foundation; either version 2 of the License, or
0015 //  (at your option) any later version.
0016 //
0017 //  The algorithms that underlie FastJet have required considerable
0018 //  development. They are described in the original FastJet paper,
0019 //  hep-ph/0512210 and in the manual, arXiv:1111.6097. If you use
0020 //  FastJet as part of work towards a scientific publication, please
0021 //  quote the version you use and include a citation to the manual and
0022 //  optionally also to hep-ph/0512210.
0023 //
0024 //  FastJet is distributed in the hope that it will be useful,
0025 //  but WITHOUT ANY WARRANTY; without even the implied warranty of
0026 //  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
0027 //  GNU General Public License for more details.
0028 //
0029 //  You should have received a copy of the GNU General Public License
0030 //  along with FastJet. If not, see <http://www.gnu.org/licenses/>.
0031 //----------------------------------------------------------------------
0032 //FJENDHEADER
0033
0034 //#include "fastjet/PseudoJet.hh"
0035 #include "fastjet/internal/MinHeap.hh"
0036 #include "fastjet/ClusterSequence.hh"
0037
0038 FASTJET_BEGIN_NAMESPACE      // defined in fastjet/internal/base.hh
0039
0040 /// Rounding errors in the Lazy strategies may cause the following
0041 /// problem: when browsing tiles in the vicinity of the particles
0042 /// being clustered in order to decide which of these tiles may
0043 /// contain particles that need to be updated (because theit NN is one
0044 /// of the particles that are currently clustered), we discard tiles
0045 /// that are deemed "too far from the cell" by the "max_NN_dist"
0046 /// criterion. Because of rounding error, this condition can sometimes
0047 /// miss cases where an update is needed.
0048 ///
0049 /// An example of this happens if a particle '1' is, say, at the lower
0050 /// edge of the rapidity of a given tile, with a particle '2' in the
0051 /// tile directly on its left at the same rapidity. Assume also that
0052 /// max_NN_dist in 2's tile corresponds to the distance between 2 and
0053 /// teh tile of 1. If 2 is 1's NN then in case 2 gets clustered, 1's
0054 /// NN needs to be updated. However, rounding errors in the
0055 /// calculation of the distance between 1 and 2 may result is
0056 /// something slightly larger than the max_NN_dist in 2's tile.
0057 ///
0058 /// This situation corresponds to the bug reported by Jochen Olt on
0059 /// February 12 2015 [see issue-tracker/2015-02-infinite-loop],
0060 /// causing an infinite loop.
0061 ///
0062 /// To prevent this, the simplest solution is, when looking at tiles
0063 /// to browse for updateds, to add a margin of security close to the
0064 /// edges of the cell, i.e. instead of updating only tiles for which
0065 /// distance<=max_NN_dist, we will update tiles for which
0066 /// distance<=max_NN_dist+tile_edge_security_margin.
0067 ///
0068 /// Note that this does not need to be done when computing nearest
0069 /// neighbours [rounding errors are tolerated there] but it is
0070 /// critical when tracking points that have to be updated.
0071 const double tile_edge_security_margin=1.0e-7;
0072
0073 /// structure analogous to BriefJet, but with the extra information
0074 /// needed for dealing with tiles
0075 class TiledJet {
0076 public:
0077   double     eta, phi, kt2, NN_dist;
0078   TiledJet * NN, *previous, * next;
0079   int        _jets_index, tile_index;
0080   bool _minheap_update_needed;
0081
0082   // indicate whether jets need to have their minheap entries
0083   // updated).
0084   inline void label_minheap_update_needed() {_minheap_update_needed = true;}
0085   inline void label_minheap_update_done()   {_minheap_update_needed = false;}
0086   inline bool minheap_update_needed() const {return _minheap_update_needed;}
0087 };
0088
0089 const int n_tile_neighbours = 9;
0090
0091 class Tile {
0092 public:
0093   typedef double (Tile::*DistToTileFn)(const TiledJet*) const;
0094   typedef std::pair<Tile *, DistToTileFn> TileFnPair;
0095   /// pointers to neighbouring tiles, including self
0096   TileFnPair begin_tiles[n_tile_neighbours];
0097   /// neighbouring tiles, excluding self
0098   TileFnPair *  surrounding_tiles;
0099   /// half of neighbouring tiles, no self
0100   TileFnPair *  RH_tiles;
0101   /// just beyond end of tiles
0102   TileFnPair *  end_tiles;
0103   /// start of list of BriefJets contained in this tile
0104   TiledJet * head;
0105   /// sometimes useful to be able to tag a tile
0106   bool     tagged;
0107   /// true for tiles where the delta phi calculation needs
0108   /// potentially to account for periodicity in phi
0109   bool     use_periodic_delta_phi;
0110   /// for all particles in the tile, this stores the largest of the
0111   /// (squared) nearest-neighbour distances.
0112   double max_NN_dist;
0113   double eta_min, eta_max, phi_min, phi_max;
0114
0115   double distance_to_centre(const TiledJet *) const {return 0;}
0116   double distance_to_left(const TiledJet * jet) const {
0117     double deta = jet->eta - eta_min;
0118     return deta*deta;
0119   }
0120   double distance_to_right(const TiledJet * jet) const {
0121     double deta = jet->eta - eta_max;
0122     return deta*deta;
0123   }
0124   double distance_to_bottom(const TiledJet * jet) const {
0125     double dphi = jet->phi - phi_min;
0126     return dphi*dphi;
0127   }
0128   double distance_to_top(const TiledJet * jet) const {
0129     double dphi = jet->phi - phi_max;
0130     return dphi*dphi;
0131   }
0132
0133   double distance_to_left_top(const TiledJet * jet) const {
0134     double deta = jet->eta - eta_min;
0135     double dphi = jet->phi - phi_max;
0136     return deta*deta + dphi*dphi;
0137   }
0138   double distance_to_left_bottom(const TiledJet * jet) const {
0139     double deta = jet->eta - eta_min;
0140     double dphi = jet->phi - phi_min;
0141     return deta*deta + dphi*dphi;
0142   }
0143   double distance_to_right_top(const TiledJet * jet) const {
0144     double deta = jet->eta - eta_max;
0145     double dphi = jet->phi - phi_max;
0146     return deta*deta + dphi*dphi;
0147   }
0148   double distance_to_right_bottom(const TiledJet * jet) const {
0149     double deta = jet->eta - eta_max;
0150     double dphi = jet->phi - phi_min;
0151     return deta*deta + dphi*dphi;
0152   }
0153
0154
0155 };
0156
0157 //----------------------------------------------------------------------
0158 class LazyTiling9Alt {
0159 public:
0160   LazyTiling9Alt(ClusterSequence & cs);
0161
0162   void run();
0163
0164   //void get_next_clustering(int & jetA_index, int & jetB_index, double & dij);
0165
0166
0167 protected:
0168   ClusterSequence & _cs;
0169   const std::vector<PseudoJet> & _jets;
0170   std::vector<Tile> _tiles;
0171
0172
0173   double _Rparam, _R2, _invR2;
0174   double _tiles_eta_min, _tiles_eta_max;
0175   double _tile_size_eta, _tile_size_phi;
0176   double _tile_half_size_eta, _tile_half_size_phi;
0177   int    _n_tiles_phi,_tiles_ieta_min,_tiles_ieta_max;
0178
0179   std::vector<TiledJet *> _jets_for_minheap;
0180
0181   //MinHeap _minheap;
0182
0183   void _initialise_tiles();
0184
0185   // reasonably robust return of tile index given ieta and iphi, in particular
0186   // it works even if iphi is negative
0187   inline int _tile_index (int ieta, int iphi) const {
0188     // note that (-1)%n = -1 so that we have to add _n_tiles_phi
0189     // before performing modulo operation
0190     return (ieta-_tiles_ieta_min)*_n_tiles_phi
0191                   + (iphi+_n_tiles_phi) % _n_tiles_phi;
0192   }
0193
0194   void  _bj_remove_from_tiles(TiledJet * const jet);
0195
0196   /// returns the tile index given the eta and phi values of a jet
0197   int _tile_index(const double eta, const double phi) const;
0198
0199   // sets up information regarding the tiling of the given jet
0200   void _tj_set_jetinfo(TiledJet * const jet, const int _jets_index);
0201
0202   void _print_tiles(TiledJet * briefjets ) const;
0203   void _add_neighbours_to_tile_union(const int tile_index,
0204          std::vector<int> & tile_union, int & n_near_tiles) const;
0205   void _add_untagged_neighbours_to_tile_union(const int tile_index,
0206          std::vector<int> & tile_union, int & n_near_tiles);
0207   void _add_untagged_neighbours_to_tile_union_using_max_info(const TiledJet * const jet,
0208          std::vector<int> & tile_union, int & n_near_tiles);
0209   //double _distance_to_tile(const TiledJet * bj, const Tile *) const;
0210   void _update_jetX_jetI_NN(TiledJet * jetX, TiledJet * jetI, std::vector<TiledJet *> & jets_for_minheap);
0211
0212   void _set_NN(TiledJet * jetI, std::vector<TiledJet *> & jets_for_minheap);
0213
0214   // return the diJ (multiplied by _R2) for this jet assuming its NN
0215   // info is correct
0216   template <class J> double _bj_diJ(const J * const jet) const {
0217     double kt2 = jet->kt2;
0218     if (jet->NN != NULL) {if (jet->NN->kt2 < kt2) {kt2 = jet->NN->kt2;}}
0219     return jet->NN_dist * kt2;
0220   }
0221
0222
0223   //----------------------------------------------------------------------
0224   template <class J> inline void _bj_set_jetinfo(
0225                             J * const jetA, const int _jets_index) const {
0226     jetA->eta  = _jets[_jets_index].rap();
0227     jetA->phi  = _jets[_jets_index].phi_02pi();
0228     jetA->kt2  = _cs.jet_scale_for_algorithm(_jets[_jets_index]);
0229     jetA->_jets_index = _jets_index;
0230     // initialise NN info as well
0231     jetA->NN_dist = _R2;
0232     jetA->NN      = NULL;
0233   }
0234
0235
0236   //----------------------------------------------------------------------
0237   template <class J> inline double _bj_dist(
0238                 const J * const jetA, const J * const jetB) const {
0239     double dphi = std::abs(jetA->phi - jetB->phi);
0240     double deta = (jetA->eta - jetB->eta);
0241     if (dphi > pi) {dphi = twopi - dphi;}
0242     return dphi*dphi + deta*deta;
0243   }
0244
0245
0246   //----------------------------------------------------------------------
0247   template <class J> inline double _bj_dist_not_periodic(
0248                 const J * const jetA, const J * const jetB) const {
0249     double dphi = jetA->phi - jetB->phi;
0250     double deta = (jetA->eta - jetB->eta);
0251     return dphi*dphi + deta*deta;
0252   }
0253
0254 };
0255
0256
0257 FASTJET_END_NAMESPACE
0258
0259 #endif // __FASTJET_LAZYTILING9ALT_HH__