FastJet  3.1.0-beta.1
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Groups Pages
LazyTiling9Alt.hh
1 #ifndef __FASTJET_LAZYTILING9ALT_HH__
2 #define __FASTJET_LAZYTILING9ALT_HH__
3 
4 //FJSTARTHEADER
5 // $Id: LazyTiling9Alt.hh 3477 2014-07-29 14:34:39Z salam $
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 
38 FASTJET_BEGIN_NAMESPACE // defined in fastjet/internal/base.hh
39 
40 /// structure analogous to BriefJet, but with the extra information
41 /// needed for dealing with tiles
42 class TiledJet {
43 public:
44  double eta, phi, kt2, NN_dist;
45  TiledJet * NN, *previous, * next;
46  int _jets_index, tile_index;
47  bool _minheap_update_needed;
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 const int n_tile_neighbours = 9;
57 
58 class Tile {
59 public:
60  typedef double (Tile::*DistToTileFn)(const TiledJet*) const;
61  typedef std::pair<Tile *, DistToTileFn> TileFnPair;
62  /// pointers to neighbouring tiles, including self
63  TileFnPair begin_tiles[n_tile_neighbours];
64  /// neighbouring tiles, excluding self
65  TileFnPair * surrounding_tiles;
66  /// half of neighbouring tiles, no self
67  TileFnPair * RH_tiles;
68  /// just beyond end of tiles
69  TileFnPair * end_tiles;
70  /// start of list of BriefJets contained in this tile
71  TiledJet * head;
72  /// sometimes useful to be able to tag a tile
73  bool tagged;
74  /// true for tiles where the delta phi calculation needs
75  /// potentially to account for periodicity in phi
76  bool use_periodic_delta_phi;
77  /// for all particles in the tile, this stores the largest of the
78  /// (squared) nearest-neighbour distances.
79  double max_NN_dist;
80  double eta_min, eta_max, phi_min, phi_max;
81 
82  double distance_to_centre(const TiledJet *) const {return 0;}
83  double distance_to_left(const TiledJet * jet) const {
84  double deta = jet->eta - eta_min;
85  return deta*deta;
86  }
87  double distance_to_right(const TiledJet * jet) const {
88  double deta = jet->eta - eta_max;
89  return deta*deta;
90  }
91  double distance_to_bottom(const TiledJet * jet) const {
92  double dphi = jet->phi - phi_min;
93  return dphi*dphi;
94  }
95  double distance_to_top(const TiledJet * jet) const {
96  double dphi = jet->phi - phi_max;
97  return dphi*dphi;
98  }
99 
100  double distance_to_left_top(const TiledJet * jet) const {
101  double deta = jet->eta - eta_min;
102  double dphi = jet->phi - phi_max;
103  return deta*deta + dphi*dphi;
104  }
105  double distance_to_left_bottom(const TiledJet * jet) const {
106  double deta = jet->eta - eta_min;
107  double dphi = jet->phi - phi_min;
108  return deta*deta + dphi*dphi;
109  }
110  double distance_to_right_top(const TiledJet * jet) const {
111  double deta = jet->eta - eta_max;
112  double dphi = jet->phi - phi_max;
113  return deta*deta + dphi*dphi;
114  }
115  double distance_to_right_bottom(const TiledJet * jet) const {
116  double deta = jet->eta - eta_max;
117  double dphi = jet->phi - phi_min;
118  return deta*deta + dphi*dphi;
119  }
120 
121 
122 };
123 
124 //----------------------------------------------------------------------
125 class LazyTiling9Alt {
126 public:
127  LazyTiling9Alt(ClusterSequence & cs);
128 
129  void run();
130 
131  //void get_next_clustering(int & jetA_index, int & jetB_index, double & dij);
132 
133 
134 protected:
135  ClusterSequence & _cs;
136  const std::vector<PseudoJet> & _jets;
137  std::vector<Tile> _tiles;
138 
139 
140  double _Rparam, _R2, _invR2;
141  double _tiles_eta_min, _tiles_eta_max;
142  double _tile_size_eta, _tile_size_phi;
143  double _tile_half_size_eta, _tile_half_size_phi;
144  int _n_tiles_phi,_tiles_ieta_min,_tiles_ieta_max;
145 
146  std::vector<TiledJet *> _jets_for_minheap;
147 
148  //MinHeap _minheap;
149 
150  void _initialise_tiles();
151 
152  // reasonably robust return of tile index given ieta and iphi, in particular
153  // it works even if iphi is negative
154  inline int _tile_index (int ieta, int iphi) const {
155  // note that (-1)%n = -1 so that we have to add _n_tiles_phi
156  // before performing modulo operation
157  return (ieta-_tiles_ieta_min)*_n_tiles_phi
158  + (iphi+_n_tiles_phi) % _n_tiles_phi;
159  }
160 
161  void _bj_remove_from_tiles(TiledJet * const jet);
162 
163  /// returns the tile index given the eta and phi values of a jet
164  int _tile_index(const double eta, const double phi) const;
165 
166  // sets up information regarding the tiling of the given jet
167  void _tj_set_jetinfo(TiledJet * const jet, const int _jets_index);
168 
169  void _print_tiles(TiledJet * briefjets ) const;
170  void _add_neighbours_to_tile_union(const int tile_index,
171  std::vector<int> & tile_union, int & n_near_tiles) const;
172  void _add_untagged_neighbours_to_tile_union(const int tile_index,
173  std::vector<int> & tile_union, int & n_near_tiles);
174  void _add_untagged_neighbours_to_tile_union_using_max_info(const TiledJet * const jet,
175  std::vector<int> & tile_union, int & n_near_tiles);
176  //double _distance_to_tile(const TiledJet * bj, const Tile *) const;
177  void _update_jetX_jetI_NN(TiledJet * jetX, TiledJet * jetI, std::vector<TiledJet *> & jets_for_minheap);
178 
179  void _set_NN(TiledJet * jetI, std::vector<TiledJet *> & jets_for_minheap);
180 
181  // return the diJ (multiplied by _R2) for this jet assuming its NN
182  // info is correct
183  template <class J> double _bj_diJ(const J * const jet) const {
184  double kt2 = jet->kt2;
185  if (jet->NN != NULL) {if (jet->NN->kt2 < kt2) {kt2 = jet->NN->kt2;}}
186  return jet->NN_dist * kt2;
187  }
188 
189 
190  //----------------------------------------------------------------------
191  template <class J> inline void _bj_set_jetinfo(
192  J * const jetA, const int _jets_index) const {
193  jetA->eta = _jets[_jets_index].rap();
194  jetA->phi = _jets[_jets_index].phi_02pi();
195  jetA->kt2 = _cs.jet_scale_for_algorithm(_jets[_jets_index]);
196  jetA->_jets_index = _jets_index;
197  // initialise NN info as well
198  jetA->NN_dist = _R2;
199  jetA->NN = NULL;
200  }
201 
202 
203  //----------------------------------------------------------------------
204  template <class J> inline double _bj_dist(
205  const J * const jetA, const J * const jetB) const {
206  double dphi = std::abs(jetA->phi - jetB->phi);
207  double deta = (jetA->eta - jetB->eta);
208  if (dphi > pi) {dphi = twopi - dphi;}
209  return dphi*dphi + deta*deta;
210  }
211 
212 
213  //----------------------------------------------------------------------
214  template <class J> inline double _bj_dist_not_periodic(
215  const J * const jetA, const J * const jetB) const {
216  double dphi = jetA->phi - jetB->phi;
217  double deta = (jetA->eta - jetB->eta);
218  return dphi*dphi + deta*deta;
219  }
220 
221 };
222 
223 
224 FASTJET_END_NAMESPACE
225 
226 #endif // __FASTJET_LAZYTILING9ALT_HH__