fastjet/internal/LazyTiling9SeparateGhosts.hh

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