Fork me on GitHub

source: svn/trunk/Utilities/Fastjet/src/ClusterSequence_TiledN2.cc@ 99

Last change on this file since 99 was 11, checked in by severine ovyn, 16 years ago

Fastjet added; CDFCones directory has been changed

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