fastjet/internal/LazyTiling25.hh

0001 #ifndef __FASTJET_LAZYTILING25_HH__
0002 #define __FASTJET_LAZYTILING25_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 #include "fastjet/internal/LazyTiling9.hh"
0041
0042
0043
0044 FASTJET_BEGIN_NAMESPACE      // defined in fastjet/internal/base.hh
0045
0046 typedef Tile2Base<25> Tile25;
0047
0048
0049 // class Tile25 {
0050 // public:
0051 //   /// pointers to neighbouring tiles, including self
0052 //   Tile25 *   begin_tiles[25];
0053 //   /// neighbouring tiles, excluding self
0054 //   Tile25 **  surrounding_tiles;
0055 //   /// half of neighbouring tiles, no self
0056 //   Tile25 **  RH_tiles;
0057 //   /// just beyond end of tiles
0058 //   Tile25 **  end_tiles;
0059 //   /// start of list of BriefJets contained in this tile
0060 //   TiledJet * head;
0061 //   /// sometimes useful to be able to tag a tile
0062 //   bool     tagged;
0063 //   /// for all particles in the tile, this stores the largest of the
0064 //   /// (squared) nearest-neighbour distances.
0065 //   double max_NN_dist;
0066 //   double eta_centre, phi_centre;
0067 //
0068 //   bool is_near_zero_phi(double tile_size_phi) const {
0069 //     return phi_centre < 2*tile_size_phi || (twopi-phi_centre) < 2*tile_size_phi;
0070 //   }
0071 // };
0072
0073
0074 //----------------------------------------------------------------------
0075 class LazyTiling25 {
0076 public:
0077   LazyTiling25(ClusterSequence & cs);
0078
0079   void run();
0080
0081 protected:
0082   ClusterSequence & _cs;
0083   const std::vector<PseudoJet> & _jets;
0084   std::vector<Tile25> _tiles;
0085
0086 #ifdef INSTRUMENT2
0087   int _ncall; // GPS tmp
0088   int _ncall_dtt; // GPS tmp
0089 #endif // INSTRUMENT2
0090
0091   double _Rparam, _R2, _invR2;
0092   double _tiles_eta_min, _tiles_eta_max;
0093   double _tile_size_eta, _tile_size_phi;
0094   double _tile_half_size_eta, _tile_half_size_phi;
0095   int    _n_tiles_phi,_tiles_ieta_min,_tiles_ieta_max;
0096
0097   std::vector<TiledJet *> _jets_for_minheap;
0098
0099   //MinHeap _minheap;
0100
0101   void _initialise_tiles();
0102
0103   // reasonably robust return of tile index given ieta and iphi, in particular
0104   // it works even if iphi is negative
0105   inline int _tile_index (int ieta, int iphi) const {
0106     // note that (-1)%n = -1 so that we have to add _n_tiles_phi
0107     // before performing modulo operation
0108     return (ieta-_tiles_ieta_min)*_n_tiles_phi
0109                   + (iphi+_n_tiles_phi) % _n_tiles_phi;
0110   }
0111
0112   void  _bj_remove_from_tiles(TiledJet * const jet);
0113
0114   /// returns the tile index given the eta and phi values of a jet
0115   int _tile_index(const double eta, const double phi) const;
0116
0117   // sets up information regarding the tiling of the given jet
0118   void _tj_set_jetinfo(TiledJet * const jet, const int _jets_index);
0119
0120   void _print_tiles(TiledJet * briefjets ) const;
0121   void _add_neighbours_to_tile_union(const int tile_index,
0122          std::vector<int> & tile_union, int & n_near_tiles) const;
0123   void _add_untagged_neighbours_to_tile_union(const int tile_index,
0124          std::vector<int> & tile_union, int & n_near_tiles);
0125   void _add_untagged_neighbours_to_tile_union_using_max_info(const TiledJet * const jet,
0126          std::vector<int> & tile_union, int & n_near_tiles);
0127   double _distance_to_tile(const TiledJet * bj, const Tile25 *)
0128 #ifdef INSTRUMENT2
0129     ;
0130 #else
0131     const;
0132 #endif
0133   void _update_jetX_jetI_NN(TiledJet * jetX, TiledJet * jetI, std::vector<TiledJet *> & jets_for_minheap);
0134
0135   void _set_NN(TiledJet * jetI, std::vector<TiledJet *> & jets_for_minheap);
0136
0137   // return the diJ (multiplied by _R2) for this jet assuming its NN
0138   // info is correct
0139   template <class J> double _bj_diJ(const J * const jet) const {
0140     double kt2 = jet->kt2;
0141     if (jet->NN != NULL) {if (jet->NN->kt2 < kt2) {kt2 = jet->NN->kt2;}}
0142     return jet->NN_dist * kt2;
0143   }
0144
0145
0146   //----------------------------------------------------------------------
0147   template <class J> inline void _bj_set_jetinfo(
0148                             J * const jetA, const int _jets_index) const {
0149     jetA->eta  = _jets[_jets_index].rap();
0150     jetA->phi  = _jets[_jets_index].phi_02pi();
0151     jetA->kt2  = _cs.jet_scale_for_algorithm(_jets[_jets_index]);
0152     jetA->_jets_index = _jets_index;
0153     // initialise NN info as well
0154     jetA->NN_dist = _R2;
0155     jetA->NN      = NULL;
0156   }
0157
0158
0159   //----------------------------------------------------------------------
0160   template <class J> inline double _bj_dist(
0161                 const J * const jetA, const J * const jetB)
0162 #ifdef INSTRUMENT2
0163     {
0164     _ncall++; // GPS tmp
0165 #else
0166     const {
0167 #endif
0168     double dphi = std::abs(jetA->phi - jetB->phi);
0169     double deta = (jetA->eta - jetB->eta);
0170     if (dphi > pi) {dphi = twopi - dphi;}
0171     return dphi*dphi + deta*deta;
0172   }
0173
0174
0175   //----------------------------------------------------------------------
0176   template <class J> inline double _bj_dist_not_periodic(
0177                 const J * const jetA, const J * const jetB)
0178 #ifdef INSTRUMENT2
0179     {
0180     _ncall++; // GPS tmp
0181 #else
0182     const {
0183 #endif
0184     //_ncall++; // GPS tmp
0185     double dphi = jetA->phi - jetB->phi;
0186     double deta = (jetA->eta - jetB->eta);
0187     return dphi*dphi + deta*deta;
0188   }
0189
0190 };
0191
0192
0193 FASTJET_END_NAMESPACE
0194
0195 #endif // __FASTJET_LAZYTILING25_HH__