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