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