1 | //STARTHEADER |
---|
2 | // $Id: ClusterSequence_TiledN2.cc 859 2012-11-28 01:49:23Z pavel $ |
---|
3 | // |
---|
4 | // Copyright (c) 2005-2011, Matteo Cacciari, Gavin P. Salam and Gregory Soyez |
---|
5 | // |
---|
6 | //---------------------------------------------------------------------- |
---|
7 | // This file is part of FastJet. |
---|
8 | // |
---|
9 | // FastJet is free software; you can redistribute it and/or modify |
---|
10 | // it under the terms of the GNU General Public License as published by |
---|
11 | // the Free Software Foundation; either version 2 of the License, or |
---|
12 | // (at your option) any later version. |
---|
13 | // |
---|
14 | // The algorithms that underlie FastJet have required considerable |
---|
15 | // development and are described in hep-ph/0512210. If you use |
---|
16 | // FastJet as part of work towards a scientific publication, please |
---|
17 | // include a citation to the FastJet paper. |
---|
18 | // |
---|
19 | // FastJet is distributed in the hope that it will be useful, |
---|
20 | // but WITHOUT ANY WARRANTY; without even the implied warranty of |
---|
21 | // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the |
---|
22 | // GNU General Public License for more details. |
---|
23 | // |
---|
24 | // You should have received a copy of the GNU General Public License |
---|
25 | // along with FastJet. If not, see <http://www.gnu.org/licenses/>. |
---|
26 | //---------------------------------------------------------------------- |
---|
27 | //ENDHEADER |
---|
28 | |
---|
29 | |
---|
30 | // The plain N^2 part of the ClusterSequence class -- separated out |
---|
31 | // from the rest of the class implementation so as to speed up |
---|
32 | // compilation of this particular part while it is under test. |
---|
33 | |
---|
34 | #include<iostream> |
---|
35 | #include<vector> |
---|
36 | #include<cmath> |
---|
37 | #include<algorithm> |
---|
38 | #include "fastjet/PseudoJet.hh" |
---|
39 | #include "fastjet/ClusterSequence.hh" |
---|
40 | #include "fastjet/internal/MinHeap.hh" |
---|
41 | |
---|
42 | FASTJET_BEGIN_NAMESPACE // defined in fastjet/internal/base.hh |
---|
43 | |
---|
44 | using namespace std; |
---|
45 | |
---|
46 | |
---|
47 | //---------------------------------------------------------------------- |
---|
48 | void ClusterSequence::_bj_remove_from_tiles(TiledJet * const jet) { |
---|
49 | Tile * tile = & _tiles[jet->tile_index]; |
---|
50 | |
---|
51 | if (jet->previous == NULL) { |
---|
52 | // we are at head of the tile, so reset it. |
---|
53 | // If this was the only jet on the tile then tile->head will now be NULL |
---|
54 | tile->head = jet->next; |
---|
55 | } else { |
---|
56 | // adjust link from previous jet in this tile |
---|
57 | jet->previous->next = jet->next; |
---|
58 | } |
---|
59 | if (jet->next != NULL) { |
---|
60 | // adjust backwards-link from next jet in this tile |
---|
61 | jet->next->previous = jet->previous; |
---|
62 | } |
---|
63 | } |
---|
64 | |
---|
65 | //---------------------------------------------------------------------- |
---|
66 | /// Set up the tiles: |
---|
67 | /// - decide the range in eta |
---|
68 | /// - allocate the tiles |
---|
69 | /// - set up the cross-referencing info between tiles |
---|
70 | /// |
---|
71 | /// The neighbourhood of a tile is set up as follows |
---|
72 | /// |
---|
73 | /// LRR |
---|
74 | /// LXR |
---|
75 | /// LLR |
---|
76 | /// |
---|
77 | /// such that tiles is an array containing XLLLLRRRR with pointers |
---|
78 | /// | \ RH_tiles |
---|
79 | /// \ surrounding_tiles |
---|
80 | /// |
---|
81 | /// with appropriate precautions when close to the edge of the tiled |
---|
82 | /// region. |
---|
83 | /// |
---|
84 | void ClusterSequence::_initialise_tiles() { |
---|
85 | |
---|
86 | // first decide tile sizes (with a lower bound to avoid huge memory use with |
---|
87 | // very small R) |
---|
88 | double default_size = max(0.1,_Rparam); |
---|
89 | _tile_size_eta = default_size; |
---|
90 | // it makes no sense to go below 3 tiles in phi -- 3 tiles is |
---|
91 | // sufficient to make sure all pair-wise combinations up to pi in |
---|
92 | // phi are possible |
---|
93 | _n_tiles_phi = max(3,int(floor(twopi/default_size))); |
---|
94 | _tile_size_phi = twopi / _n_tiles_phi; // >= _Rparam and fits in 2pi |
---|
95 | |
---|
96 | // always include zero rapidity in the tiling region |
---|
97 | _tiles_eta_min = 0.0; |
---|
98 | _tiles_eta_max = 0.0; |
---|
99 | // but go no further than following |
---|
100 | const double maxrap = 7.0; |
---|
101 | |
---|
102 | // and find out how much further one should go |
---|
103 | for(unsigned int i = 0; i < _jets.size(); i++) { |
---|
104 | double eta = _jets[i].rap(); |
---|
105 | // first check if eta is in range -- to avoid taking into account |
---|
106 | // very spurious rapidities due to particles with near-zero kt. |
---|
107 | if (abs(eta) < maxrap) { |
---|
108 | if (eta < _tiles_eta_min) {_tiles_eta_min = eta;} |
---|
109 | if (eta > _tiles_eta_max) {_tiles_eta_max = eta;} |
---|
110 | } |
---|
111 | } |
---|
112 | |
---|
113 | // now adjust the values |
---|
114 | _tiles_ieta_min = int(floor(_tiles_eta_min/_tile_size_eta)); |
---|
115 | _tiles_ieta_max = int(floor( _tiles_eta_max/_tile_size_eta)); |
---|
116 | _tiles_eta_min = _tiles_ieta_min * _tile_size_eta; |
---|
117 | _tiles_eta_max = _tiles_ieta_max * _tile_size_eta; |
---|
118 | |
---|
119 | // allocate the tiles |
---|
120 | _tiles.resize((_tiles_ieta_max-_tiles_ieta_min+1)*_n_tiles_phi); |
---|
121 | |
---|
122 | // now set up the cross-referencing between tiles |
---|
123 | for (int ieta = _tiles_ieta_min; ieta <= _tiles_ieta_max; ieta++) { |
---|
124 | for (int iphi = 0; iphi < _n_tiles_phi; iphi++) { |
---|
125 | Tile * tile = & _tiles[_tile_index(ieta,iphi)]; |
---|
126 | // no jets in this tile yet |
---|
127 | tile->head = NULL; // first element of tiles points to itself |
---|
128 | tile->begin_tiles[0] = tile; |
---|
129 | Tile ** pptile = & (tile->begin_tiles[0]); |
---|
130 | pptile++; |
---|
131 | // |
---|
132 | // set up L's in column to the left of X |
---|
133 | tile->surrounding_tiles = pptile; |
---|
134 | if (ieta > _tiles_ieta_min) { |
---|
135 | // with the itile subroutine, we can safely run tiles from |
---|
136 | // idphi=-1 to idphi=+1, because it takes care of |
---|
137 | // negative and positive boundaries |
---|
138 | for (int idphi = -1; idphi <=+1; idphi++) { |
---|
139 | *pptile = & _tiles[_tile_index(ieta-1,iphi+idphi)]; |
---|
140 | pptile++; |
---|
141 | } |
---|
142 | } |
---|
143 | // now set up last L (below X) |
---|
144 | *pptile = & _tiles[_tile_index(ieta,iphi-1)]; |
---|
145 | pptile++; |
---|
146 | // set up first R (above X) |
---|
147 | tile->RH_tiles = pptile; |
---|
148 | *pptile = & _tiles[_tile_index(ieta,iphi+1)]; |
---|
149 | pptile++; |
---|
150 | // set up remaining R's, to the right of X |
---|
151 | if (ieta < _tiles_ieta_max) { |
---|
152 | for (int idphi = -1; idphi <= +1; idphi++) { |
---|
153 | *pptile = & _tiles[_tile_index(ieta+1,iphi+idphi)]; |
---|
154 | pptile++; |
---|
155 | } |
---|
156 | } |
---|
157 | // now put semaphore for end tile |
---|
158 | tile->end_tiles = pptile; |
---|
159 | // finally make sure tiles are untagged |
---|
160 | tile->tagged = false; |
---|
161 | } |
---|
162 | } |
---|
163 | |
---|
164 | } |
---|
165 | |
---|
166 | |
---|
167 | //---------------------------------------------------------------------- |
---|
168 | /// return the tile index corresponding to the given eta,phi point |
---|
169 | int ClusterSequence::_tile_index(const double & eta, const double & phi) const { |
---|
170 | int ieta, iphi; |
---|
171 | if (eta <= _tiles_eta_min) {ieta = 0;} |
---|
172 | else if (eta >= _tiles_eta_max) {ieta = _tiles_ieta_max-_tiles_ieta_min;} |
---|
173 | else { |
---|
174 | //ieta = int(floor((eta - _tiles_eta_min) / _tile_size_eta)); |
---|
175 | ieta = int(((eta - _tiles_eta_min) / _tile_size_eta)); |
---|
176 | // following needed in case of rare but nasty rounding errors |
---|
177 | if (ieta > _tiles_ieta_max-_tiles_ieta_min) { |
---|
178 | ieta = _tiles_ieta_max-_tiles_ieta_min;} |
---|
179 | } |
---|
180 | // allow for some extent of being beyond range in calculation of phi |
---|
181 | // as well |
---|
182 | //iphi = (int(floor(phi/_tile_size_phi)) + _n_tiles_phi) % _n_tiles_phi; |
---|
183 | // with just int and no floor, things run faster but beware |
---|
184 | iphi = int((phi+twopi)/_tile_size_phi) % _n_tiles_phi; |
---|
185 | return (iphi + ieta * _n_tiles_phi); |
---|
186 | } |
---|
187 | |
---|
188 | |
---|
189 | //---------------------------------------------------------------------- |
---|
190 | // overloaded version which additionally sets up information regarding the |
---|
191 | // tiling |
---|
192 | inline void ClusterSequence::_tj_set_jetinfo( TiledJet * const jet, |
---|
193 | const int _jets_index) { |
---|
194 | // first call the generic setup |
---|
195 | _bj_set_jetinfo<>(jet, _jets_index); |
---|
196 | |
---|
197 | // Then do the setup specific to the tiled case. |
---|
198 | |
---|
199 | // Find out which tile it belonds to |
---|
200 | jet->tile_index = _tile_index(jet->eta, jet->phi); |
---|
201 | |
---|
202 | // Insert it into the tile's linked list of jets |
---|
203 | Tile * tile = &_tiles[jet->tile_index]; |
---|
204 | jet->previous = NULL; |
---|
205 | jet->next = tile->head; |
---|
206 | if (jet->next != NULL) {jet->next->previous = jet;} |
---|
207 | tile->head = jet; |
---|
208 | } |
---|
209 | |
---|
210 | |
---|
211 | //---------------------------------------------------------------------- |
---|
212 | /// output the contents of the tiles |
---|
213 | void ClusterSequence::_print_tiles(TiledJet * briefjets ) const { |
---|
214 | for (vector<Tile>::const_iterator tile = _tiles.begin(); |
---|
215 | tile < _tiles.end(); tile++) { |
---|
216 | cout << "Tile " << tile - _tiles.begin()<<" = "; |
---|
217 | vector<int> list; |
---|
218 | for (TiledJet * jetI = tile->head; jetI != NULL; jetI = jetI->next) { |
---|
219 | list.push_back(jetI-briefjets); |
---|
220 | //cout <<" "<<jetI-briefjets; |
---|
221 | } |
---|
222 | sort(list.begin(),list.end()); |
---|
223 | for (unsigned int i = 0; i < list.size(); i++) {cout <<" "<<list[i];} |
---|
224 | cout <<"\n"; |
---|
225 | } |
---|
226 | } |
---|
227 | |
---|
228 | |
---|
229 | //---------------------------------------------------------------------- |
---|
230 | /// Add to the vector tile_union the tiles that are in the neighbourhood |
---|
231 | /// of the specified tile_index, including itself -- start adding |
---|
232 | /// from position n_near_tiles-1, and increase n_near_tiles as |
---|
233 | /// you go along (could have done it more C++ like with vector with reserved |
---|
234 | /// space, but fear is that it would have been slower, e.g. checking |
---|
235 | /// for end of vector at each stage to decide whether to resize it) |
---|
236 | void ClusterSequence::_add_neighbours_to_tile_union(const int tile_index, |
---|
237 | vector<int> & tile_union, int & n_near_tiles) const { |
---|
238 | for (Tile * const * near_tile = _tiles[tile_index].begin_tiles; |
---|
239 | near_tile != _tiles[tile_index].end_tiles; near_tile++){ |
---|
240 | // get the tile number |
---|
241 | tile_union[n_near_tiles] = *near_tile - & _tiles[0]; |
---|
242 | n_near_tiles++; |
---|
243 | } |
---|
244 | } |
---|
245 | |
---|
246 | |
---|
247 | //---------------------------------------------------------------------- |
---|
248 | /// Like _add_neighbours_to_tile_union, but only adds neighbours if |
---|
249 | /// their "tagged" status is false; when a neighbour is added its |
---|
250 | /// tagged status is set to true. |
---|
251 | inline void ClusterSequence::_add_untagged_neighbours_to_tile_union( |
---|
252 | const int tile_index, |
---|
253 | vector<int> & tile_union, int & n_near_tiles) { |
---|
254 | for (Tile ** near_tile = _tiles[tile_index].begin_tiles; |
---|
255 | near_tile != _tiles[tile_index].end_tiles; near_tile++){ |
---|
256 | if (! (*near_tile)->tagged) { |
---|
257 | (*near_tile)->tagged = true; |
---|
258 | // get the tile number |
---|
259 | tile_union[n_near_tiles] = *near_tile - & _tiles[0]; |
---|
260 | n_near_tiles++; |
---|
261 | } |
---|
262 | } |
---|
263 | } |
---|
264 | |
---|
265 | |
---|
266 | //---------------------------------------------------------------------- |
---|
267 | /// run a tiled clustering |
---|
268 | void ClusterSequence::_tiled_N2_cluster() { |
---|
269 | |
---|
270 | _initialise_tiles(); |
---|
271 | |
---|
272 | int n = _jets.size(); |
---|
273 | TiledJet * briefjets = new TiledJet[n]; |
---|
274 | TiledJet * jetA = briefjets, * jetB; |
---|
275 | TiledJet oldB; |
---|
276 | |
---|
277 | |
---|
278 | // will be used quite deep inside loops, but declare it here so that |
---|
279 | // memory (de)allocation gets done only once |
---|
280 | vector<int> tile_union(3*n_tile_neighbours); |
---|
281 | |
---|
282 | // initialise the basic jet info |
---|
283 | for (int i = 0; i< n; i++) { |
---|
284 | _tj_set_jetinfo(jetA, i); |
---|
285 | //cout << i<<": "<<jetA->tile_index<<"\n"; |
---|
286 | jetA++; // move on to next entry of briefjets |
---|
287 | } |
---|
288 | TiledJet * tail = jetA; // a semaphore for the end of briefjets |
---|
289 | TiledJet * head = briefjets; // a nicer way of naming start |
---|
290 | |
---|
291 | // set up the initial nearest neighbour information |
---|
292 | vector<Tile>::const_iterator tile; |
---|
293 | for (tile = _tiles.begin(); tile != _tiles.end(); tile++) { |
---|
294 | // first do it on this tile |
---|
295 | for (jetA = tile->head; jetA != NULL; jetA = jetA->next) { |
---|
296 | for (jetB = tile->head; jetB != jetA; jetB = jetB->next) { |
---|
297 | double dist = _bj_dist(jetA,jetB); |
---|
298 | if (dist < jetA->NN_dist) {jetA->NN_dist = dist; jetA->NN = jetB;} |
---|
299 | if (dist < jetB->NN_dist) {jetB->NN_dist = dist; jetB->NN = jetA;} |
---|
300 | } |
---|
301 | } |
---|
302 | // then do it for RH tiles |
---|
303 | for (Tile ** RTile = tile->RH_tiles; RTile != tile->end_tiles; RTile++) { |
---|
304 | for (jetA = tile->head; jetA != NULL; jetA = jetA->next) { |
---|
305 | for (jetB = (*RTile)->head; jetB != NULL; jetB = jetB->next) { |
---|
306 | double dist = _bj_dist(jetA,jetB); |
---|
307 | if (dist < jetA->NN_dist) {jetA->NN_dist = dist; jetA->NN = jetB;} |
---|
308 | if (dist < jetB->NN_dist) {jetB->NN_dist = dist; jetB->NN = jetA;} |
---|
309 | } |
---|
310 | } |
---|
311 | } |
---|
312 | } |
---|
313 | |
---|
314 | // now create the diJ (where J is i's NN) table -- remember that |
---|
315 | // we differ from standard normalisation here by a factor of R2 |
---|
316 | double * diJ = new double[n]; |
---|
317 | jetA = head; |
---|
318 | for (int i = 0; i < n; i++) { |
---|
319 | diJ[i] = _bj_diJ(jetA); |
---|
320 | jetA++; // have jetA follow i |
---|
321 | } |
---|
322 | |
---|
323 | // now run the recombination loop |
---|
324 | int history_location = n-1; |
---|
325 | while (tail != head) { |
---|
326 | |
---|
327 | // find the minimum of the diJ on this round |
---|
328 | double diJ_min = diJ[0]; |
---|
329 | int diJ_min_jet = 0; |
---|
330 | for (int i = 1; i < n; i++) { |
---|
331 | if (diJ[i] < diJ_min) {diJ_min_jet = i; diJ_min = diJ[i];} |
---|
332 | } |
---|
333 | |
---|
334 | // do the recombination between A and B |
---|
335 | history_location++; |
---|
336 | jetA = & briefjets[diJ_min_jet]; |
---|
337 | jetB = jetA->NN; |
---|
338 | // put the normalisation back in |
---|
339 | diJ_min *= _invR2; |
---|
340 | |
---|
341 | //if (n == 19) {cout << "Hello "<<jetA-head<<" "<<jetB-head<<"\n";} |
---|
342 | |
---|
343 | //cout <<" WILL RECOMBINE "<< jetA-briefjets<<" "<<jetB-briefjets<<"\n"; |
---|
344 | |
---|
345 | if (jetB != NULL) { |
---|
346 | // jet-jet recombination |
---|
347 | // If necessary relabel A & B to ensure jetB < jetA, that way if |
---|
348 | // the larger of them == newtail then that ends up being jetA and |
---|
349 | // the new jet that is added as jetB is inserted in a position that |
---|
350 | // has a future! |
---|
351 | if (jetA < jetB) {std::swap(jetA,jetB);} |
---|
352 | |
---|
353 | int nn; // new jet index |
---|
354 | _do_ij_recombination_step(jetA->_jets_index, jetB->_jets_index, diJ_min, nn); |
---|
355 | |
---|
356 | //OBS// get the two history indices |
---|
357 | //OBSint hist_a = _jets[jetA->_jets_index].cluster_hist_index(); |
---|
358 | //OBSint hist_b = _jets[jetB->_jets_index].cluster_hist_index(); |
---|
359 | //OBS// create the recombined jet |
---|
360 | //OBS_jets.push_back(_jets[jetA->_jets_index] + _jets[jetB->_jets_index]); |
---|
361 | //OBSint nn = _jets.size() - 1; |
---|
362 | //OBS_jets[nn].set_cluster_hist_index(history_location); |
---|
363 | //OBS// update history |
---|
364 | //OBS//cout <<n-1<<" "<<jetA-head<<" "<<jetB-head<<"; "; |
---|
365 | //OBS_add_step_to_history(history_location, |
---|
366 | //OBS min(hist_a,hist_b),max(hist_a,hist_b), |
---|
367 | //OBS nn, diJ_min); |
---|
368 | |
---|
369 | // what was jetB will now become the new jet |
---|
370 | _bj_remove_from_tiles(jetA); |
---|
371 | oldB = * jetB; // take a copy because we will need it... |
---|
372 | _bj_remove_from_tiles(jetB); |
---|
373 | _tj_set_jetinfo(jetB, nn); // also registers the jet in the tiling |
---|
374 | } else { |
---|
375 | // jet-beam recombination |
---|
376 | _do_iB_recombination_step(jetA->_jets_index, diJ_min); |
---|
377 | |
---|
378 | //OBS// get the hist_index |
---|
379 | //OBSint hist_a = _jets[jetA->_jets_index].cluster_hist_index(); |
---|
380 | //OBS//cout <<n-1<<" "<<jetA-head<<" "<<-1<<"; "; |
---|
381 | //OBS_add_step_to_history(history_location,hist_a,BeamJet,Invalid,diJ_min); |
---|
382 | _bj_remove_from_tiles(jetA); |
---|
383 | } |
---|
384 | |
---|
385 | // first establish the set of tiles over which we are going to |
---|
386 | // have to run searches for updated and new nearest-neighbours |
---|
387 | int n_near_tiles = 0; |
---|
388 | _add_neighbours_to_tile_union(jetA->tile_index, tile_union, n_near_tiles); |
---|
389 | if (jetB != NULL) { |
---|
390 | bool sort_it = false; |
---|
391 | if (jetB->tile_index != jetA->tile_index) { |
---|
392 | sort_it = true; |
---|
393 | _add_neighbours_to_tile_union(jetB->tile_index,tile_union,n_near_tiles); |
---|
394 | } |
---|
395 | if (oldB.tile_index != jetA->tile_index && |
---|
396 | oldB.tile_index != jetB->tile_index) { |
---|
397 | sort_it = true; |
---|
398 | _add_neighbours_to_tile_union(oldB.tile_index,tile_union,n_near_tiles); |
---|
399 | } |
---|
400 | |
---|
401 | if (sort_it) { |
---|
402 | // sort the tiles before then compressing the list |
---|
403 | sort(tile_union.begin(), tile_union.begin()+n_near_tiles); |
---|
404 | // and now condense the list |
---|
405 | int nnn = 1; |
---|
406 | for (int i = 1; i < n_near_tiles; i++) { |
---|
407 | if (tile_union[i] != tile_union[nnn-1]) { |
---|
408 | tile_union[nnn] = tile_union[i]; |
---|
409 | nnn++; |
---|
410 | } |
---|
411 | } |
---|
412 | n_near_tiles = nnn; |
---|
413 | } |
---|
414 | } |
---|
415 | |
---|
416 | // now update our nearest neighbour info and diJ table |
---|
417 | // first reduce size of table |
---|
418 | tail--; n--; |
---|
419 | if (jetA == tail) { |
---|
420 | // there is nothing to be done |
---|
421 | } else { |
---|
422 | // Copy last jet contents and diJ info into position of jetA |
---|
423 | *jetA = *tail; |
---|
424 | diJ[jetA - head] = diJ[tail-head]; |
---|
425 | // IN the tiling fix pointers to tail and turn them into |
---|
426 | // pointers to jetA (from predecessors, successors and the tile |
---|
427 | // head if need be) |
---|
428 | if (jetA->previous == NULL) { |
---|
429 | _tiles[jetA->tile_index].head = jetA; |
---|
430 | } else { |
---|
431 | jetA->previous->next = jetA; |
---|
432 | } |
---|
433 | if (jetA->next != NULL) {jetA->next->previous = jetA;} |
---|
434 | } |
---|
435 | |
---|
436 | // Initialise jetB's NN distance as well as updating it for |
---|
437 | // other particles. |
---|
438 | for (int itile = 0; itile < n_near_tiles; itile++) { |
---|
439 | Tile * tile_ptr = &_tiles[tile_union[itile]]; |
---|
440 | for (TiledJet * jetI = tile_ptr->head; jetI != NULL; jetI = jetI->next) { |
---|
441 | // see if jetI had jetA or jetB as a NN -- if so recalculate the NN |
---|
442 | if (jetI->NN == jetA || (jetI->NN == jetB && jetB != NULL)) { |
---|
443 | jetI->NN_dist = _R2; |
---|
444 | jetI->NN = NULL; |
---|
445 | // now go over tiles that are neighbours of I (include own tile) |
---|
446 | for (Tile ** near_tile = tile_ptr->begin_tiles; |
---|
447 | near_tile != tile_ptr->end_tiles; near_tile++) { |
---|
448 | // and then over the contents of that tile |
---|
449 | for (TiledJet * jetJ = (*near_tile)->head; |
---|
450 | jetJ != NULL; jetJ = jetJ->next) { |
---|
451 | double dist = _bj_dist(jetI,jetJ); |
---|
452 | if (dist < jetI->NN_dist && jetJ != jetI) { |
---|
453 | jetI->NN_dist = dist; jetI->NN = jetJ; |
---|
454 | } |
---|
455 | } |
---|
456 | } |
---|
457 | diJ[jetI-head] = _bj_diJ(jetI); // update diJ |
---|
458 | } |
---|
459 | // check whether new jetB is closer than jetI's current NN and |
---|
460 | // if need to update things |
---|
461 | if (jetB != NULL) { |
---|
462 | double dist = _bj_dist(jetI,jetB); |
---|
463 | if (dist < jetI->NN_dist) { |
---|
464 | if (jetI != jetB) { |
---|
465 | jetI->NN_dist = dist; |
---|
466 | jetI->NN = jetB; |
---|
467 | diJ[jetI-head] = _bj_diJ(jetI); // update diJ... |
---|
468 | } |
---|
469 | } |
---|
470 | if (dist < jetB->NN_dist) { |
---|
471 | if (jetI != jetB) { |
---|
472 | jetB->NN_dist = dist; |
---|
473 | jetB->NN = jetI;} |
---|
474 | } |
---|
475 | } |
---|
476 | } |
---|
477 | } |
---|
478 | |
---|
479 | |
---|
480 | if (jetB != NULL) {diJ[jetB-head] = _bj_diJ(jetB);} |
---|
481 | //cout << n<<" "<<briefjets[95].NN-briefjets<<" "<<briefjets[95].NN_dist <<"\n"; |
---|
482 | |
---|
483 | // remember to update pointers to tail |
---|
484 | for (Tile ** near_tile = _tiles[tail->tile_index].begin_tiles; |
---|
485 | near_tile!= _tiles[tail->tile_index].end_tiles; near_tile++){ |
---|
486 | // and then the contents of that tile |
---|
487 | for (TiledJet * jetJ = (*near_tile)->head; |
---|
488 | jetJ != NULL; jetJ = jetJ->next) { |
---|
489 | if (jetJ->NN == tail) {jetJ->NN = jetA;} |
---|
490 | } |
---|
491 | } |
---|
492 | |
---|
493 | //for (int i = 0; i < n; i++) { |
---|
494 | // if (briefjets[i].NN-briefjets >= n && briefjets[i].NN != NULL) {cout <<"YOU MUST BE CRAZY for n ="<<n<<", i = "<<i<<", NN = "<<briefjets[i].NN-briefjets<<"\n";} |
---|
495 | //} |
---|
496 | |
---|
497 | |
---|
498 | if (jetB != NULL) {diJ[jetB-head] = _bj_diJ(jetB);} |
---|
499 | //cout << briefjets[95].NN-briefjets<<" "<<briefjets[95].NN_dist <<"\n"; |
---|
500 | |
---|
501 | } |
---|
502 | |
---|
503 | // final cleaning up; |
---|
504 | delete[] diJ; |
---|
505 | delete[] briefjets; |
---|
506 | } |
---|
507 | |
---|
508 | |
---|
509 | //---------------------------------------------------------------------- |
---|
510 | /// run a tiled clustering |
---|
511 | void ClusterSequence::_faster_tiled_N2_cluster() { |
---|
512 | |
---|
513 | _initialise_tiles(); |
---|
514 | |
---|
515 | int n = _jets.size(); |
---|
516 | TiledJet * briefjets = new TiledJet[n]; |
---|
517 | TiledJet * jetA = briefjets, * jetB; |
---|
518 | TiledJet oldB; |
---|
519 | |
---|
520 | |
---|
521 | // will be used quite deep inside loops, but declare it here so that |
---|
522 | // memory (de)allocation gets done only once |
---|
523 | vector<int> tile_union(3*n_tile_neighbours); |
---|
524 | |
---|
525 | // initialise the basic jet info |
---|
526 | for (int i = 0; i< n; i++) { |
---|
527 | _tj_set_jetinfo(jetA, i); |
---|
528 | //cout << i<<": "<<jetA->tile_index<<"\n"; |
---|
529 | jetA++; // move on to next entry of briefjets |
---|
530 | } |
---|
531 | TiledJet * head = briefjets; // a nicer way of naming start |
---|
532 | |
---|
533 | // set up the initial nearest neighbour information |
---|
534 | vector<Tile>::const_iterator tile; |
---|
535 | for (tile = _tiles.begin(); tile != _tiles.end(); tile++) { |
---|
536 | // first do it on this tile |
---|
537 | for (jetA = tile->head; jetA != NULL; jetA = jetA->next) { |
---|
538 | for (jetB = tile->head; jetB != jetA; jetB = jetB->next) { |
---|
539 | double dist = _bj_dist(jetA,jetB); |
---|
540 | if (dist < jetA->NN_dist) {jetA->NN_dist = dist; jetA->NN = jetB;} |
---|
541 | if (dist < jetB->NN_dist) {jetB->NN_dist = dist; jetB->NN = jetA;} |
---|
542 | } |
---|
543 | } |
---|
544 | // then do it for RH tiles |
---|
545 | for (Tile ** RTile = tile->RH_tiles; RTile != tile->end_tiles; RTile++) { |
---|
546 | for (jetA = tile->head; jetA != NULL; jetA = jetA->next) { |
---|
547 | for (jetB = (*RTile)->head; jetB != NULL; jetB = jetB->next) { |
---|
548 | double dist = _bj_dist(jetA,jetB); |
---|
549 | if (dist < jetA->NN_dist) {jetA->NN_dist = dist; jetA->NN = jetB;} |
---|
550 | if (dist < jetB->NN_dist) {jetB->NN_dist = dist; jetB->NN = jetA;} |
---|
551 | } |
---|
552 | } |
---|
553 | } |
---|
554 | // no need to do it for LH tiles, since they are implicitly done |
---|
555 | // when we set NN for both jetA and jetB on the RH tiles. |
---|
556 | } |
---|
557 | |
---|
558 | |
---|
559 | // now create the diJ (where J is i's NN) table -- remember that |
---|
560 | // we differ from standard normalisation here by a factor of R2 |
---|
561 | // (corrected for at the end). |
---|
562 | struct diJ_plus_link { |
---|
563 | double diJ; // the distance |
---|
564 | TiledJet * jet; // the jet (i) for which we've found this distance |
---|
565 | // (whose NN will the J). |
---|
566 | }; |
---|
567 | diJ_plus_link * diJ = new diJ_plus_link[n]; |
---|
568 | jetA = head; |
---|
569 | for (int i = 0; i < n; i++) { |
---|
570 | diJ[i].diJ = _bj_diJ(jetA); // kt distance * R^2 |
---|
571 | diJ[i].jet = jetA; // our compact diJ table will not be in |
---|
572 | jetA->diJ_posn = i; // one-to-one corresp. with non-compact jets, |
---|
573 | // so set up bi-directional correspondence here. |
---|
574 | jetA++; // have jetA follow i |
---|
575 | } |
---|
576 | |
---|
577 | // now run the recombination loop |
---|
578 | int history_location = n-1; |
---|
579 | while (n > 0) { |
---|
580 | |
---|
581 | // find the minimum of the diJ on this round |
---|
582 | diJ_plus_link * best, *stop; // pointers a bit faster than indices |
---|
583 | // could use best to keep track of diJ min, but it turns out to be |
---|
584 | // marginally faster to have a separate variable (avoids n |
---|
585 | // dereferences at the expense of n/2 assignments). |
---|
586 | double diJ_min = diJ[0].diJ; // initialise the best one here. |
---|
587 | best = diJ; // and here |
---|
588 | stop = diJ+n; |
---|
589 | for (diJ_plus_link * here = diJ+1; here != stop; here++) { |
---|
590 | if (here->diJ < diJ_min) {best = here; diJ_min = here->diJ;} |
---|
591 | } |
---|
592 | |
---|
593 | // do the recombination between A and B |
---|
594 | history_location++; |
---|
595 | jetA = best->jet; |
---|
596 | jetB = jetA->NN; |
---|
597 | // put the normalisation back in |
---|
598 | diJ_min *= _invR2; |
---|
599 | |
---|
600 | if (jetB != NULL) { |
---|
601 | // jet-jet recombination |
---|
602 | // If necessary relabel A & B to ensure jetB < jetA, that way if |
---|
603 | // the larger of them == newtail then that ends up being jetA and |
---|
604 | // the new jet that is added as jetB is inserted in a position that |
---|
605 | // has a future! |
---|
606 | if (jetA < jetB) {std::swap(jetA,jetB);} |
---|
607 | |
---|
608 | int nn; // new jet index |
---|
609 | _do_ij_recombination_step(jetA->_jets_index, jetB->_jets_index, diJ_min, nn); |
---|
610 | |
---|
611 | //OBS// get the two history indices |
---|
612 | //OBSint ihstry_a = _jets[jetA->_jets_index].cluster_hist_index(); |
---|
613 | //OBSint ihstry_b = _jets[jetB->_jets_index].cluster_hist_index(); |
---|
614 | //OBS// create the recombined jet |
---|
615 | //OBS_jets.push_back(_jets[jetA->_jets_index] + _jets[jetB->_jets_index]); |
---|
616 | //OBSint nn = _jets.size() - 1; |
---|
617 | //OBS_jets[nn].set_cluster_hist_index(history_location); |
---|
618 | //OBS// update history |
---|
619 | //OBS//cout <<n-1<<" "<<jetA-head<<" "<<jetB-head<<"; "; |
---|
620 | //OBS_add_step_to_history(history_location, |
---|
621 | //OBS min(ihstry_a,ihstry_b),max(ihstry_a,ihstry_b), |
---|
622 | //OBS nn, diJ_min); |
---|
623 | // what was jetB will now become the new jet |
---|
624 | _bj_remove_from_tiles(jetA); |
---|
625 | oldB = * jetB; // take a copy because we will need it... |
---|
626 | _bj_remove_from_tiles(jetB); |
---|
627 | _tj_set_jetinfo(jetB, nn); // cause jetB to become _jets[nn] |
---|
628 | // (also registers the jet in the tiling) |
---|
629 | } else { |
---|
630 | // jet-beam recombination |
---|
631 | // get the hist_index |
---|
632 | _do_iB_recombination_step(jetA->_jets_index, diJ_min); |
---|
633 | //OBSint ihstry_a = _jets[jetA->_jets_index].cluster_hist_index(); |
---|
634 | //OBS//cout <<n-1<<" "<<jetA-head<<" "<<-1<<"; "; |
---|
635 | //OBS_add_step_to_history(history_location,ihstry_a,BeamJet,Invalid,diJ_min); |
---|
636 | _bj_remove_from_tiles(jetA); |
---|
637 | } |
---|
638 | |
---|
639 | // first establish the set of tiles over which we are going to |
---|
640 | // have to run searches for updated and new nearest-neighbours -- |
---|
641 | // basically a combination of vicinity of the tiles of the two old |
---|
642 | // and one new jet. |
---|
643 | int n_near_tiles = 0; |
---|
644 | _add_untagged_neighbours_to_tile_union(jetA->tile_index, |
---|
645 | tile_union, n_near_tiles); |
---|
646 | if (jetB != NULL) { |
---|
647 | if (jetB->tile_index != jetA->tile_index) { |
---|
648 | _add_untagged_neighbours_to_tile_union(jetB->tile_index, |
---|
649 | tile_union,n_near_tiles); |
---|
650 | } |
---|
651 | if (oldB.tile_index != jetA->tile_index && |
---|
652 | oldB.tile_index != jetB->tile_index) { |
---|
653 | _add_untagged_neighbours_to_tile_union(oldB.tile_index, |
---|
654 | tile_union,n_near_tiles); |
---|
655 | } |
---|
656 | } |
---|
657 | |
---|
658 | // now update our nearest neighbour info and diJ table |
---|
659 | // first reduce size of table |
---|
660 | n--; |
---|
661 | // then compactify the diJ by taking the last of the diJ and copying |
---|
662 | // it to the position occupied by the diJ for jetA |
---|
663 | diJ[n].jet->diJ_posn = jetA->diJ_posn; |
---|
664 | diJ[jetA->diJ_posn] = diJ[n]; |
---|
665 | |
---|
666 | // Initialise jetB's NN distance as well as updating it for |
---|
667 | // other particles. |
---|
668 | // Run over all tiles in our union |
---|
669 | for (int itile = 0; itile < n_near_tiles; itile++) { |
---|
670 | Tile * tile_ptr = &_tiles[tile_union[itile]]; |
---|
671 | tile_ptr->tagged = false; // reset tag, since we're done with unions |
---|
672 | // run over all jets in the current tile |
---|
673 | for (TiledJet * jetI = tile_ptr->head; jetI != NULL; jetI = jetI->next) { |
---|
674 | // see if jetI had jetA or jetB as a NN -- if so recalculate the NN |
---|
675 | if (jetI->NN == jetA || (jetI->NN == jetB && jetB != NULL)) { |
---|
676 | jetI->NN_dist = _R2; |
---|
677 | jetI->NN = NULL; |
---|
678 | // now go over tiles that are neighbours of I (include own tile) |
---|
679 | for (Tile ** near_tile = tile_ptr->begin_tiles; |
---|
680 | near_tile != tile_ptr->end_tiles; near_tile++) { |
---|
681 | // and then over the contents of that tile |
---|
682 | for (TiledJet * jetJ = (*near_tile)->head; |
---|
683 | jetJ != NULL; jetJ = jetJ->next) { |
---|
684 | double dist = _bj_dist(jetI,jetJ); |
---|
685 | if (dist < jetI->NN_dist && jetJ != jetI) { |
---|
686 | jetI->NN_dist = dist; jetI->NN = jetJ; |
---|
687 | } |
---|
688 | } |
---|
689 | } |
---|
690 | diJ[jetI->diJ_posn].diJ = _bj_diJ(jetI); // update diJ kt-dist |
---|
691 | } |
---|
692 | // check whether new jetB is closer than jetI's current NN and |
---|
693 | // if jetI is closer than jetB's current (evolving) nearest |
---|
694 | // neighbour. Where relevant update things |
---|
695 | if (jetB != NULL) { |
---|
696 | double dist = _bj_dist(jetI,jetB); |
---|
697 | if (dist < jetI->NN_dist) { |
---|
698 | if (jetI != jetB) { |
---|
699 | jetI->NN_dist = dist; |
---|
700 | jetI->NN = jetB; |
---|
701 | diJ[jetI->diJ_posn].diJ = _bj_diJ(jetI); // update diJ... |
---|
702 | } |
---|
703 | } |
---|
704 | if (dist < jetB->NN_dist) { |
---|
705 | if (jetI != jetB) { |
---|
706 | jetB->NN_dist = dist; |
---|
707 | jetB->NN = jetI;} |
---|
708 | } |
---|
709 | } |
---|
710 | } |
---|
711 | } |
---|
712 | |
---|
713 | // finally, register the updated kt distance for B |
---|
714 | if (jetB != NULL) {diJ[jetB->diJ_posn].diJ = _bj_diJ(jetB);} |
---|
715 | |
---|
716 | } |
---|
717 | |
---|
718 | // final cleaning up; |
---|
719 | delete[] diJ; |
---|
720 | delete[] briefjets; |
---|
721 | } |
---|
722 | |
---|
723 | |
---|
724 | |
---|
725 | //---------------------------------------------------------------------- |
---|
726 | /// run a tiled clustering, with our minheap for keeping track of the |
---|
727 | /// smallest dij |
---|
728 | void ClusterSequence::_minheap_faster_tiled_N2_cluster() { |
---|
729 | |
---|
730 | _initialise_tiles(); |
---|
731 | |
---|
732 | int n = _jets.size(); |
---|
733 | TiledJet * briefjets = new TiledJet[n]; |
---|
734 | TiledJet * jetA = briefjets, * jetB; |
---|
735 | TiledJet oldB; |
---|
736 | |
---|
737 | |
---|
738 | // will be used quite deep inside loops, but declare it here so that |
---|
739 | // memory (de)allocation gets done only once |
---|
740 | vector<int> tile_union(3*n_tile_neighbours); |
---|
741 | |
---|
742 | // initialise the basic jet info |
---|
743 | for (int i = 0; i< n; i++) { |
---|
744 | _tj_set_jetinfo(jetA, i); |
---|
745 | //cout << i<<": "<<jetA->tile_index<<"\n"; |
---|
746 | jetA++; // move on to next entry of briefjets |
---|
747 | } |
---|
748 | TiledJet * head = briefjets; // a nicer way of naming start |
---|
749 | |
---|
750 | // set up the initial nearest neighbour information |
---|
751 | vector<Tile>::const_iterator tile; |
---|
752 | for (tile = _tiles.begin(); tile != _tiles.end(); tile++) { |
---|
753 | // first do it on this tile |
---|
754 | for (jetA = tile->head; jetA != NULL; jetA = jetA->next) { |
---|
755 | for (jetB = tile->head; jetB != jetA; jetB = jetB->next) { |
---|
756 | double dist = _bj_dist(jetA,jetB); |
---|
757 | if (dist < jetA->NN_dist) {jetA->NN_dist = dist; jetA->NN = jetB;} |
---|
758 | if (dist < jetB->NN_dist) {jetB->NN_dist = dist; jetB->NN = jetA;} |
---|
759 | } |
---|
760 | } |
---|
761 | // then do it for RH tiles |
---|
762 | for (Tile ** RTile = tile->RH_tiles; RTile != tile->end_tiles; RTile++) { |
---|
763 | for (jetA = tile->head; jetA != NULL; jetA = jetA->next) { |
---|
764 | for (jetB = (*RTile)->head; jetB != NULL; jetB = jetB->next) { |
---|
765 | double dist = _bj_dist(jetA,jetB); |
---|
766 | if (dist < jetA->NN_dist) {jetA->NN_dist = dist; jetA->NN = jetB;} |
---|
767 | if (dist < jetB->NN_dist) {jetB->NN_dist = dist; jetB->NN = jetA;} |
---|
768 | } |
---|
769 | } |
---|
770 | } |
---|
771 | // no need to do it for LH tiles, since they are implicitly done |
---|
772 | // when we set NN for both jetA and jetB on the RH tiles. |
---|
773 | } |
---|
774 | |
---|
775 | |
---|
776 | //// now create the diJ (where J is i's NN) table -- remember that |
---|
777 | //// we differ from standard normalisation here by a factor of R2 |
---|
778 | //// (corrected for at the end). |
---|
779 | //struct diJ_plus_link { |
---|
780 | // double diJ; // the distance |
---|
781 | // TiledJet * jet; // the jet (i) for which we've found this distance |
---|
782 | // // (whose NN will the J). |
---|
783 | //}; |
---|
784 | //diJ_plus_link * diJ = new diJ_plus_link[n]; |
---|
785 | //jetA = head; |
---|
786 | //for (int i = 0; i < n; i++) { |
---|
787 | // diJ[i].diJ = _bj_diJ(jetA); // kt distance * R^2 |
---|
788 | // diJ[i].jet = jetA; // our compact diJ table will not be in |
---|
789 | // jetA->diJ_posn = i; // one-to-one corresp. with non-compact jets, |
---|
790 | // // so set up bi-directional correspondence here. |
---|
791 | // jetA++; // have jetA follow i |
---|
792 | //} |
---|
793 | |
---|
794 | vector<double> diJs(n); |
---|
795 | for (int i = 0; i < n; i++) { |
---|
796 | diJs[i] = _bj_diJ(&briefjets[i]); |
---|
797 | briefjets[i].label_minheap_update_done(); |
---|
798 | } |
---|
799 | MinHeap minheap(diJs); |
---|
800 | // have a stack telling us which jets we'll have to update on the heap |
---|
801 | vector<TiledJet *> jets_for_minheap; |
---|
802 | jets_for_minheap.reserve(n); |
---|
803 | |
---|
804 | // now run the recombination loop |
---|
805 | int history_location = n-1; |
---|
806 | while (n > 0) { |
---|
807 | |
---|
808 | double diJ_min = minheap.minval() *_invR2; |
---|
809 | jetA = head + minheap.minloc(); |
---|
810 | |
---|
811 | // do the recombination between A and B |
---|
812 | history_location++; |
---|
813 | jetB = jetA->NN; |
---|
814 | |
---|
815 | if (jetB != NULL) { |
---|
816 | // jet-jet recombination |
---|
817 | // If necessary relabel A & B to ensure jetB < jetA, that way if |
---|
818 | // the larger of them == newtail then that ends up being jetA and |
---|
819 | // the new jet that is added as jetB is inserted in a position that |
---|
820 | // has a future! |
---|
821 | if (jetA < jetB) {std::swap(jetA,jetB);} |
---|
822 | |
---|
823 | int nn; // new jet index |
---|
824 | _do_ij_recombination_step(jetA->_jets_index, jetB->_jets_index, diJ_min, nn); |
---|
825 | |
---|
826 | // what was jetB will now become the new jet |
---|
827 | _bj_remove_from_tiles(jetA); |
---|
828 | oldB = * jetB; // take a copy because we will need it... |
---|
829 | _bj_remove_from_tiles(jetB); |
---|
830 | _tj_set_jetinfo(jetB, nn); // cause jetB to become _jets[nn] |
---|
831 | // (also registers the jet in the tiling) |
---|
832 | } else { |
---|
833 | // jet-beam recombination |
---|
834 | // get the hist_index |
---|
835 | _do_iB_recombination_step(jetA->_jets_index, diJ_min); |
---|
836 | _bj_remove_from_tiles(jetA); |
---|
837 | } |
---|
838 | |
---|
839 | // remove the minheap entry for jetA |
---|
840 | minheap.remove(jetA-head); |
---|
841 | |
---|
842 | // first establish the set of tiles over which we are going to |
---|
843 | // have to run searches for updated and new nearest-neighbours -- |
---|
844 | // basically a combination of vicinity of the tiles of the two old |
---|
845 | // and one new jet. |
---|
846 | int n_near_tiles = 0; |
---|
847 | _add_untagged_neighbours_to_tile_union(jetA->tile_index, |
---|
848 | tile_union, n_near_tiles); |
---|
849 | if (jetB != NULL) { |
---|
850 | if (jetB->tile_index != jetA->tile_index) { |
---|
851 | _add_untagged_neighbours_to_tile_union(jetB->tile_index, |
---|
852 | tile_union,n_near_tiles); |
---|
853 | } |
---|
854 | if (oldB.tile_index != jetA->tile_index && |
---|
855 | oldB.tile_index != jetB->tile_index) { |
---|
856 | // GS: the line below generates a warning that oldB.tile_index |
---|
857 | // may be used uninitialised. However, to reach this point, we |
---|
858 | // ned jetB != NULL (see test a few lines above) and is jetB |
---|
859 | // !=NULL, one would have gone through "oldB = *jetB before |
---|
860 | // (see piece of code ~20 line above), so the index is |
---|
861 | // initialised. We do not do anything to avoid the warning to |
---|
862 | // avoid any potential speed impact. |
---|
863 | _add_untagged_neighbours_to_tile_union(oldB.tile_index, |
---|
864 | tile_union,n_near_tiles); |
---|
865 | } |
---|
866 | // indicate that we'll have to update jetB in the minheap |
---|
867 | jetB->label_minheap_update_needed(); |
---|
868 | jets_for_minheap.push_back(jetB); |
---|
869 | } |
---|
870 | |
---|
871 | |
---|
872 | // Initialise jetB's NN distance as well as updating it for |
---|
873 | // other particles. |
---|
874 | // Run over all tiles in our union |
---|
875 | for (int itile = 0; itile < n_near_tiles; itile++) { |
---|
876 | Tile * tile_ptr = &_tiles[tile_union[itile]]; |
---|
877 | tile_ptr->tagged = false; // reset tag, since we're done with unions |
---|
878 | // run over all jets in the current tile |
---|
879 | for (TiledJet * jetI = tile_ptr->head; jetI != NULL; jetI = jetI->next) { |
---|
880 | // see if jetI had jetA or jetB as a NN -- if so recalculate the NN |
---|
881 | if (jetI->NN == jetA || (jetI->NN == jetB && jetB != NULL)) { |
---|
882 | jetI->NN_dist = _R2; |
---|
883 | jetI->NN = NULL; |
---|
884 | // label jetI as needing heap action... |
---|
885 | if (!jetI->minheap_update_needed()) { |
---|
886 | jetI->label_minheap_update_needed(); |
---|
887 | jets_for_minheap.push_back(jetI);} |
---|
888 | // now go over tiles that are neighbours of I (include own tile) |
---|
889 | for (Tile ** near_tile = tile_ptr->begin_tiles; |
---|
890 | near_tile != tile_ptr->end_tiles; near_tile++) { |
---|
891 | // and then over the contents of that tile |
---|
892 | for (TiledJet * jetJ = (*near_tile)->head; |
---|
893 | jetJ != NULL; jetJ = jetJ->next) { |
---|
894 | double dist = _bj_dist(jetI,jetJ); |
---|
895 | if (dist < jetI->NN_dist && jetJ != jetI) { |
---|
896 | jetI->NN_dist = dist; jetI->NN = jetJ; |
---|
897 | } |
---|
898 | } |
---|
899 | } |
---|
900 | } |
---|
901 | // check whether new jetB is closer than jetI's current NN and |
---|
902 | // if jetI is closer than jetB's current (evolving) nearest |
---|
903 | // neighbour. Where relevant update things |
---|
904 | if (jetB != NULL) { |
---|
905 | double dist = _bj_dist(jetI,jetB); |
---|
906 | if (dist < jetI->NN_dist) { |
---|
907 | if (jetI != jetB) { |
---|
908 | jetI->NN_dist = dist; |
---|
909 | jetI->NN = jetB; |
---|
910 | // label jetI as needing heap action... |
---|
911 | if (!jetI->minheap_update_needed()) { |
---|
912 | jetI->label_minheap_update_needed(); |
---|
913 | jets_for_minheap.push_back(jetI);} |
---|
914 | } |
---|
915 | } |
---|
916 | if (dist < jetB->NN_dist) { |
---|
917 | if (jetI != jetB) { |
---|
918 | jetB->NN_dist = dist; |
---|
919 | jetB->NN = jetI;} |
---|
920 | } |
---|
921 | } |
---|
922 | } |
---|
923 | } |
---|
924 | |
---|
925 | // deal with jets whose minheap entry needs updating |
---|
926 | while (jets_for_minheap.size() > 0) { |
---|
927 | TiledJet * jetI = jets_for_minheap.back(); |
---|
928 | jets_for_minheap.pop_back(); |
---|
929 | minheap.update(jetI-head, _bj_diJ(jetI)); |
---|
930 | jetI->label_minheap_update_done(); |
---|
931 | } |
---|
932 | n--; |
---|
933 | } |
---|
934 | |
---|
935 | // final cleaning up; |
---|
936 | delete[] briefjets; |
---|
937 | } |
---|
938 | |
---|
939 | |
---|
940 | FASTJET_END_NAMESPACE |
---|
941 | |
---|