Fork me on GitHub

source: git/external/fastjet/internal/LazyTiling9Alt.hh@ a0f5d71

3.4.3pre03
Last change on this file since a0f5d71 was b7b836a, checked in by Pavel Demin <pavel-demin@…>, 6 years ago

update FastJet library to 3.3.1 and FastJet Contrib library to 1.036

  • Property mode set to 100644
File size: 9.6 KB
RevLine 
[35cdc46]1#ifndef __FASTJET_LAZYTILING9ALT_HH__
2#define __FASTJET_LAZYTILING9ALT_HH__
3
4//FJSTARTHEADER
[b7b836a]5// $Id: LazyTiling9Alt.hh 4354 2018-04-22 07:12:37Z salam $
[35cdc46]6//
[b7b836a]7// Copyright (c) 2005-2018, Matteo Cacciari, Gavin P. Salam and Gregory Soyez
[35cdc46]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
38FASTJET_BEGIN_NAMESPACE // defined in fastjet/internal/base.hh
39
[10e33bc]40/// Rounding errors in the Lazy strategies may cause the following
41/// problem: when browsing tiles in the vicinity of the particles
42/// being clustered in order to decide which of these tiles may
43/// contain particles that need to be updated (because theit NN is one
44/// of the particles that are currently clustered), we discard tiles
45/// that are deemed "too far from the cell" by the "max_NN_dist"
46/// criterion. Because of rounding error, this condition can sometimes
47/// miss cases where an update is needed.
48///
49/// An example of this happens if a particle '1' is, say, at the lower
50/// edge of the rapidity of a given tile, with a particle '2' in the
51/// tile directly on its left at the same rapidity. Assume also that
52/// max_NN_dist in 2's tile corresponds to the distance between 2 and
53/// teh tile of 1. If 2 is 1's NN then in case 2 gets clustered, 1's
54/// NN needs to be updated. However, rounding errors in the
55/// calculation of the distance between 1 and 2 may result is
56/// something slightly larger than the max_NN_dist in 2's tile.
57///
58/// This situation corresponds to the bug reported by Jochen Olt on
59/// February 12 2015 [see issue-tracker/2015-02-infinite-loop],
60/// causing an infinite loop.
61///
62/// To prevent this, the simplest solution is, when looking at tiles
63/// to browse for updateds, to add a margin of security close to the
64/// edges of the cell, i.e. instead of updating only tiles for which
65/// distance<=max_NN_dist, we will update tiles for which
66/// distance<=max_NN_dist+tile_edge_security_margin.
67///
68/// Note that this does not need to be done when computing nearest
69/// neighbours [rounding errors are tolerated there] but it is
70/// critical when tracking points that have to be updated.
71const double tile_edge_security_margin=1.0e-7;
72
[35cdc46]73/// structure analogous to BriefJet, but with the extra information
74/// needed for dealing with tiles
75class TiledJet {
76public:
77 double eta, phi, kt2, NN_dist;
78 TiledJet * NN, *previous, * next;
79 int _jets_index, tile_index;
80 bool _minheap_update_needed;
81
82 // indicate whether jets need to have their minheap entries
83 // updated).
84 inline void label_minheap_update_needed() {_minheap_update_needed = true;}
85 inline void label_minheap_update_done() {_minheap_update_needed = false;}
86 inline bool minheap_update_needed() const {return _minheap_update_needed;}
87};
88
89const int n_tile_neighbours = 9;
90
91class Tile {
92public:
93 typedef double (Tile::*DistToTileFn)(const TiledJet*) const;
94 typedef std::pair<Tile *, DistToTileFn> TileFnPair;
95 /// pointers to neighbouring tiles, including self
96 TileFnPair begin_tiles[n_tile_neighbours];
97 /// neighbouring tiles, excluding self
98 TileFnPair * surrounding_tiles;
99 /// half of neighbouring tiles, no self
100 TileFnPair * RH_tiles;
101 /// just beyond end of tiles
102 TileFnPair * end_tiles;
103 /// start of list of BriefJets contained in this tile
104 TiledJet * head;
105 /// sometimes useful to be able to tag a tile
106 bool tagged;
107 /// true for tiles where the delta phi calculation needs
108 /// potentially to account for periodicity in phi
109 bool use_periodic_delta_phi;
110 /// for all particles in the tile, this stores the largest of the
111 /// (squared) nearest-neighbour distances.
112 double max_NN_dist;
113 double eta_min, eta_max, phi_min, phi_max;
114
115 double distance_to_centre(const TiledJet *) const {return 0;}
116 double distance_to_left(const TiledJet * jet) const {
117 double deta = jet->eta - eta_min;
118 return deta*deta;
119 }
120 double distance_to_right(const TiledJet * jet) const {
121 double deta = jet->eta - eta_max;
122 return deta*deta;
123 }
124 double distance_to_bottom(const TiledJet * jet) const {
125 double dphi = jet->phi - phi_min;
126 return dphi*dphi;
127 }
128 double distance_to_top(const TiledJet * jet) const {
129 double dphi = jet->phi - phi_max;
130 return dphi*dphi;
131 }
132
133 double distance_to_left_top(const TiledJet * jet) const {
134 double deta = jet->eta - eta_min;
135 double dphi = jet->phi - phi_max;
136 return deta*deta + dphi*dphi;
137 }
138 double distance_to_left_bottom(const TiledJet * jet) const {
139 double deta = jet->eta - eta_min;
140 double dphi = jet->phi - phi_min;
141 return deta*deta + dphi*dphi;
142 }
143 double distance_to_right_top(const TiledJet * jet) const {
144 double deta = jet->eta - eta_max;
145 double dphi = jet->phi - phi_max;
146 return deta*deta + dphi*dphi;
147 }
148 double distance_to_right_bottom(const TiledJet * jet) const {
149 double deta = jet->eta - eta_max;
150 double dphi = jet->phi - phi_min;
151 return deta*deta + dphi*dphi;
152 }
153
154
155};
156
157//----------------------------------------------------------------------
158class LazyTiling9Alt {
159public:
160 LazyTiling9Alt(ClusterSequence & cs);
161
162 void run();
163
164 //void get_next_clustering(int & jetA_index, int & jetB_index, double & dij);
165
166
167protected:
168 ClusterSequence & _cs;
169 const std::vector<PseudoJet> & _jets;
170 std::vector<Tile> _tiles;
171
172
173 double _Rparam, _R2, _invR2;
174 double _tiles_eta_min, _tiles_eta_max;
175 double _tile_size_eta, _tile_size_phi;
176 double _tile_half_size_eta, _tile_half_size_phi;
177 int _n_tiles_phi,_tiles_ieta_min,_tiles_ieta_max;
178
179 std::vector<TiledJet *> _jets_for_minheap;
180
181 //MinHeap _minheap;
182
183 void _initialise_tiles();
184
185 // reasonably robust return of tile index given ieta and iphi, in particular
186 // it works even if iphi is negative
187 inline int _tile_index (int ieta, int iphi) const {
188 // note that (-1)%n = -1 so that we have to add _n_tiles_phi
189 // before performing modulo operation
190 return (ieta-_tiles_ieta_min)*_n_tiles_phi
191 + (iphi+_n_tiles_phi) % _n_tiles_phi;
192 }
193
194 void _bj_remove_from_tiles(TiledJet * const jet);
195
196 /// returns the tile index given the eta and phi values of a jet
197 int _tile_index(const double eta, const double phi) const;
198
199 // sets up information regarding the tiling of the given jet
200 void _tj_set_jetinfo(TiledJet * const jet, const int _jets_index);
201
202 void _print_tiles(TiledJet * briefjets ) const;
203 void _add_neighbours_to_tile_union(const int tile_index,
204 std::vector<int> & tile_union, int & n_near_tiles) const;
205 void _add_untagged_neighbours_to_tile_union(const int tile_index,
206 std::vector<int> & tile_union, int & n_near_tiles);
207 void _add_untagged_neighbours_to_tile_union_using_max_info(const TiledJet * const jet,
208 std::vector<int> & tile_union, int & n_near_tiles);
209 //double _distance_to_tile(const TiledJet * bj, const Tile *) const;
210 void _update_jetX_jetI_NN(TiledJet * jetX, TiledJet * jetI, std::vector<TiledJet *> & jets_for_minheap);
211
212 void _set_NN(TiledJet * jetI, std::vector<TiledJet *> & jets_for_minheap);
213
214 // return the diJ (multiplied by _R2) for this jet assuming its NN
215 // info is correct
216 template <class J> double _bj_diJ(const J * const jet) const {
217 double kt2 = jet->kt2;
218 if (jet->NN != NULL) {if (jet->NN->kt2 < kt2) {kt2 = jet->NN->kt2;}}
219 return jet->NN_dist * kt2;
220 }
221
222
223 //----------------------------------------------------------------------
224 template <class J> inline void _bj_set_jetinfo(
225 J * const jetA, const int _jets_index) const {
226 jetA->eta = _jets[_jets_index].rap();
227 jetA->phi = _jets[_jets_index].phi_02pi();
228 jetA->kt2 = _cs.jet_scale_for_algorithm(_jets[_jets_index]);
229 jetA->_jets_index = _jets_index;
230 // initialise NN info as well
231 jetA->NN_dist = _R2;
232 jetA->NN = NULL;
233 }
234
235
236 //----------------------------------------------------------------------
237 template <class J> inline double _bj_dist(
238 const J * const jetA, const J * const jetB) const {
239 double dphi = std::abs(jetA->phi - jetB->phi);
240 double deta = (jetA->eta - jetB->eta);
241 if (dphi > pi) {dphi = twopi - dphi;}
242 return dphi*dphi + deta*deta;
243 }
244
245
246 //----------------------------------------------------------------------
247 template <class J> inline double _bj_dist_not_periodic(
248 const J * const jetA, const J * const jetB) const {
249 double dphi = jetA->phi - jetB->phi;
250 double deta = (jetA->eta - jetB->eta);
251 return dphi*dphi + deta*deta;
252 }
253
254};
255
256
257FASTJET_END_NAMESPACE
258
259#endif // __FASTJET_LAZYTILING9ALT_HH__
Note: See TracBrowser for help on using the repository browser.