| 1 | //FJSTARTHEADER
|
|---|
| 2 | // $Id: LazyTiling9SeparateGhosts.cc 3807 2015-02-20 11:16:55Z soyez $
|
|---|
| 3 | //
|
|---|
| 4 | // Copyright (c) 2005-2014, 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. They are described in the original FastJet paper,
|
|---|
| 16 | // hep-ph/0512210 and in the manual, arXiv:1111.6097. If you use
|
|---|
| 17 | // FastJet as part of work towards a scientific publication, please
|
|---|
| 18 | // quote the version you use and include a citation to the manual and
|
|---|
| 19 | // optionally also to hep-ph/0512210.
|
|---|
| 20 | //
|
|---|
| 21 | // FastJet is distributed in the hope that it will be useful,
|
|---|
| 22 | // but WITHOUT ANY WARRANTY; without even the implied warranty of
|
|---|
| 23 | // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
|
|---|
| 24 | // GNU General Public License for more details.
|
|---|
| 25 | //
|
|---|
| 26 | // You should have received a copy of the GNU General Public License
|
|---|
| 27 | // along with FastJet. If not, see <http://www.gnu.org/licenses/>.
|
|---|
| 28 | //----------------------------------------------------------------------
|
|---|
| 29 | //FJENDHEADER
|
|---|
| 30 |
|
|---|
| 31 | #include "fastjet/internal/LazyTiling9SeparateGhosts.hh"
|
|---|
| 32 | #include "fastjet/internal/TilingExtent.hh"
|
|---|
| 33 | #include <iomanip>
|
|---|
| 34 | using namespace std;
|
|---|
| 35 |
|
|---|
| 36 | FASTJET_BEGIN_NAMESPACE // defined in fastjet/internal/base.hh
|
|---|
| 37 |
|
|---|
| 38 | double LazyTiling9SeparateGhosts::ghost_pt2_threshold = 1e-100;
|
|---|
| 39 |
|
|---|
| 40 | LazyTiling9SeparateGhosts::LazyTiling9SeparateGhosts(ClusterSequence & cs) :
|
|---|
| 41 | _cs(cs), _jets(cs.jets())
|
|---|
| 42 | //, _minheap(_jets.size())
|
|---|
| 43 | {
|
|---|
| 44 | _Rparam = cs.jet_def().R();
|
|---|
| 45 | _R2 = _Rparam * _Rparam;
|
|---|
| 46 | _invR2 = 1.0 / _R2;
|
|---|
| 47 | _initialise_tiles();
|
|---|
| 48 | }
|
|---|
| 49 |
|
|---|
| 50 |
|
|---|
| 51 | //----------------------------------------------------------------------
|
|---|
| 52 | /// Set up the tiles:
|
|---|
| 53 | /// - decide the range in eta
|
|---|
| 54 | /// - allocate the tiles
|
|---|
| 55 | /// - set up the cross-referencing info between tiles
|
|---|
| 56 | ///
|
|---|
| 57 | /// The neighbourhood of a tile is set up as follows
|
|---|
| 58 | ///
|
|---|
| 59 | /// LRR
|
|---|
| 60 | /// LXR
|
|---|
| 61 | /// LLR
|
|---|
| 62 | ///
|
|---|
| 63 | /// such that tiles is an array containing XLLLLRRRR with pointers
|
|---|
| 64 | /// | \ RH_tiles
|
|---|
| 65 | /// \ surrounding_tiles
|
|---|
| 66 | ///
|
|---|
| 67 | /// with appropriate precautions when close to the edge of the tiled
|
|---|
| 68 | /// region.
|
|---|
| 69 | ///
|
|---|
| 70 | void LazyTiling9SeparateGhosts::_initialise_tiles() {
|
|---|
| 71 |
|
|---|
| 72 | // first decide tile sizes (with a lower bound to avoid huge memory use with
|
|---|
| 73 | // very small R)
|
|---|
| 74 | double default_size = max(0.1,_Rparam);
|
|---|
| 75 | _tile_size_eta = default_size;
|
|---|
| 76 | // it makes no sense to go below 3 tiles in phi -- 3 tiles is
|
|---|
| 77 | // sufficient to make sure all pair-wise combinations up to pi in
|
|---|
| 78 | // phi are possible
|
|---|
| 79 | _n_tiles_phi = max(3,int(floor(twopi/default_size)));
|
|---|
| 80 | _tile_size_phi = twopi / _n_tiles_phi; // >= _Rparam and fits in 2pi
|
|---|
| 81 |
|
|---|
| 82 | // always include zero rapidity in the tiling region
|
|---|
| 83 | _tiles_eta_min = 0.0;
|
|---|
| 84 | _tiles_eta_max = 0.0;
|
|---|
| 85 | // but go no further than following
|
|---|
| 86 | const double maxrap = 7.0;
|
|---|
| 87 |
|
|---|
| 88 | // and find out how much further one should go
|
|---|
| 89 | for(unsigned int i = 0; i < _jets.size(); i++) {
|
|---|
| 90 | double eta = _jets[i].rap();
|
|---|
| 91 | // first check if eta is in range -- to avoid taking into account
|
|---|
| 92 | // very spurious rapidities due to particles with near-zero kt.
|
|---|
| 93 | if (abs(eta) < maxrap) {
|
|---|
| 94 | if (eta < _tiles_eta_min) {_tiles_eta_min = eta;}
|
|---|
| 95 | if (eta > _tiles_eta_max) {_tiles_eta_max = eta;}
|
|---|
| 96 | }
|
|---|
| 97 | }
|
|---|
| 98 |
|
|---|
| 99 | // now adjust the values
|
|---|
| 100 | _tiles_ieta_min = int(floor(_tiles_eta_min/_tile_size_eta));
|
|---|
| 101 | _tiles_ieta_max = int(floor( _tiles_eta_max/_tile_size_eta));
|
|---|
| 102 | _tiles_eta_min = _tiles_ieta_min * _tile_size_eta;
|
|---|
| 103 | _tiles_eta_max = _tiles_ieta_max * _tile_size_eta;
|
|---|
| 104 |
|
|---|
| 105 | _tile_half_size_eta = _tile_size_eta * 0.5;
|
|---|
| 106 | _tile_half_size_phi = _tile_size_phi * 0.5;
|
|---|
| 107 |
|
|---|
| 108 | // allocate the tiles
|
|---|
| 109 | _tiles.resize((_tiles_ieta_max-_tiles_ieta_min+1)*_n_tiles_phi);
|
|---|
| 110 |
|
|---|
| 111 | // now set up the cross-referencing between tiles
|
|---|
| 112 | for (int ieta = _tiles_ieta_min; ieta <= _tiles_ieta_max; ieta++) {
|
|---|
| 113 | for (int iphi = 0; iphi < _n_tiles_phi; iphi++) {
|
|---|
| 114 | Tile3 * tile = & _tiles[_tile_index(ieta,iphi)];
|
|---|
| 115 | // no jets in this tile yet
|
|---|
| 116 | tile->head = NULL; // first element of tiles points to itself
|
|---|
| 117 | tile->ghost_head = NULL; // first element of tiles points to itself
|
|---|
| 118 | tile->begin_tiles[0] = tile;
|
|---|
| 119 | Tile3 ** pptile = & (tile->begin_tiles[0]);
|
|---|
| 120 | pptile++;
|
|---|
| 121 | //
|
|---|
| 122 | // set up L's in column to the left of X
|
|---|
| 123 | tile->surrounding_tiles = pptile;
|
|---|
| 124 | if (ieta > _tiles_ieta_min) {
|
|---|
| 125 | // with the itile subroutine, we can safely run tiles from
|
|---|
| 126 | // idphi=-1 to idphi=+1, because it takes care of
|
|---|
| 127 | // negative and positive boundaries
|
|---|
| 128 | for (int idphi = -1; idphi <=+1; idphi++) {
|
|---|
| 129 | *pptile = & _tiles[_tile_index(ieta-1,iphi+idphi)];
|
|---|
| 130 | pptile++;
|
|---|
| 131 | }
|
|---|
| 132 | }
|
|---|
| 133 | // now set up last L (below X)
|
|---|
| 134 | *pptile = & _tiles[_tile_index(ieta,iphi-1)];
|
|---|
| 135 | pptile++;
|
|---|
| 136 | // set up first R (above X)
|
|---|
| 137 | tile->RH_tiles = pptile;
|
|---|
| 138 | *pptile = & _tiles[_tile_index(ieta,iphi+1)];
|
|---|
| 139 | pptile++;
|
|---|
| 140 | // set up remaining R's, to the right of X
|
|---|
| 141 | if (ieta < _tiles_ieta_max) {
|
|---|
| 142 | for (int idphi = -1; idphi <= +1; idphi++) {
|
|---|
| 143 | *pptile = & _tiles[_tile_index(ieta+1,iphi+idphi)];
|
|---|
| 144 | pptile++;
|
|---|
| 145 | }
|
|---|
| 146 | }
|
|---|
| 147 | // now put semaphore for end tile
|
|---|
| 148 | tile->end_tiles = pptile;
|
|---|
| 149 | // finally make sure tiles are untagged
|
|---|
| 150 | tile->tagged = false;
|
|---|
| 151 | // and ensure max distance is sensibly initialised
|
|---|
| 152 | tile->max_NN_dist = 0;
|
|---|
| 153 | // and also position of centre of tile
|
|---|
| 154 | tile->eta_centre = (ieta+0.5)*_tile_size_eta;
|
|---|
| 155 | tile->phi_centre = (iphi+0.5)*_tile_size_phi;
|
|---|
| 156 | }
|
|---|
| 157 | }
|
|---|
| 158 |
|
|---|
| 159 | }
|
|---|
| 160 |
|
|---|
| 161 | //----------------------------------------------------------------------
|
|---|
| 162 | /// return the tile index corresponding to the given eta,phi point
|
|---|
| 163 | int LazyTiling9SeparateGhosts::_tile_index(const double eta, const double phi) const {
|
|---|
| 164 | int ieta, iphi;
|
|---|
| 165 | if (eta <= _tiles_eta_min) {ieta = 0;}
|
|---|
| 166 | else if (eta >= _tiles_eta_max) {ieta = _tiles_ieta_max-_tiles_ieta_min;}
|
|---|
| 167 | else {
|
|---|
| 168 | //ieta = int(floor((eta - _tiles_eta_min) / _tile_size_eta));
|
|---|
| 169 | ieta = int(((eta - _tiles_eta_min) / _tile_size_eta));
|
|---|
| 170 | // following needed in case of rare but nasty rounding errors
|
|---|
| 171 | if (ieta > _tiles_ieta_max-_tiles_ieta_min) {
|
|---|
| 172 | ieta = _tiles_ieta_max-_tiles_ieta_min;}
|
|---|
| 173 | }
|
|---|
| 174 | // allow for some extent of being beyond range in calculation of phi
|
|---|
| 175 | // as well
|
|---|
| 176 | //iphi = (int(floor(phi/_tile_size_phi)) + _n_tiles_phi) % _n_tiles_phi;
|
|---|
| 177 | // with just int and no floor, things run faster but beware
|
|---|
| 178 | iphi = int((phi+twopi)/_tile_size_phi) % _n_tiles_phi;
|
|---|
| 179 | return (iphi + ieta * _n_tiles_phi);
|
|---|
| 180 | }
|
|---|
| 181 |
|
|---|
| 182 |
|
|---|
| 183 | //----------------------------------------------------------------------
|
|---|
| 184 | // sets up information regarding the tiling of the given jet
|
|---|
| 185 | inline void LazyTiling9SeparateGhosts::_tj_set_jetinfo( TiledJet3 * const jet,
|
|---|
| 186 | const int _jets_index, bool is_ghost) {
|
|---|
| 187 | // first call the generic setup
|
|---|
| 188 | _bj_set_jetinfo<>(jet, _jets_index);
|
|---|
| 189 |
|
|---|
| 190 | // Then do the setup specific to the tiled case.
|
|---|
| 191 | jet->is_ghost = is_ghost;
|
|---|
| 192 |
|
|---|
| 193 | // Find out which tile it belonds to
|
|---|
| 194 | jet->tile_index = _tile_index(jet->eta, jet->phi);
|
|---|
| 195 |
|
|---|
| 196 | // Insert it into the tile's linked list of jets
|
|---|
| 197 | Tile3 * tile = &_tiles[jet->tile_index];
|
|---|
| 198 | jet->previous = NULL;
|
|---|
| 199 | if (is_ghost) {
|
|---|
| 200 | jet->next = tile->ghost_head;
|
|---|
| 201 | tile->ghost_head = jet;
|
|---|
| 202 | } else {
|
|---|
| 203 | jet->next = tile->head;
|
|---|
| 204 | tile->head = jet;
|
|---|
| 205 | }
|
|---|
| 206 | if (jet->next != NULL) {jet->next->previous = jet;}
|
|---|
| 207 | }
|
|---|
| 208 |
|
|---|
| 209 |
|
|---|
| 210 | //----------------------------------------------------------------------
|
|---|
| 211 | void LazyTiling9SeparateGhosts::_bj_remove_from_tiles(TiledJet3 * const jet) {
|
|---|
| 212 | Tile3 * tile = & _tiles[jet->tile_index];
|
|---|
| 213 |
|
|---|
| 214 | if (jet->previous == NULL) {
|
|---|
| 215 | // we are at head of the tile, so reset it.
|
|---|
| 216 | // If this was the only jet on the tile then tile->head will now be NULL
|
|---|
| 217 | if (jet->is_ghost) {
|
|---|
| 218 | tile->ghost_head = jet->next;
|
|---|
| 219 | } else {
|
|---|
| 220 | tile->head = jet->next;
|
|---|
| 221 | }
|
|---|
| 222 | } else {
|
|---|
| 223 | // adjust link from previous jet in this tile
|
|---|
| 224 | jet->previous->next = jet->next;
|
|---|
| 225 | }
|
|---|
| 226 | if (jet->next != NULL) {
|
|---|
| 227 | // adjust backwards-link from next jet in this tile
|
|---|
| 228 | jet->next->previous = jet->previous;
|
|---|
| 229 | }
|
|---|
| 230 | }
|
|---|
| 231 |
|
|---|
| 232 |
|
|---|
| 233 | //----------------------------------------------------------------------
|
|---|
| 234 | /// output the contents of the tiles
|
|---|
| 235 | void LazyTiling9SeparateGhosts::_print_tiles(TiledJet3 * briefjets ) const {
|
|---|
| 236 | for (vector<Tile3>::const_iterator tile = _tiles.begin();
|
|---|
| 237 | tile < _tiles.end(); tile++) {
|
|---|
| 238 | cout << "Tile " << tile - _tiles.begin()<<" = ";
|
|---|
| 239 | vector<int> list;
|
|---|
| 240 | for (TiledJet3 * jetI = tile->head; jetI != NULL; jetI = jetI->next) {
|
|---|
| 241 | list.push_back(jetI-briefjets);
|
|---|
| 242 | //cout <<" "<<jetI-briefjets;
|
|---|
| 243 | }
|
|---|
| 244 | sort(list.begin(),list.end());
|
|---|
| 245 | for (unsigned int i = 0; i < list.size(); i++) {cout <<" "<<list[i];}
|
|---|
| 246 | cout <<"\n";
|
|---|
| 247 | }
|
|---|
| 248 | }
|
|---|
| 249 |
|
|---|
| 250 |
|
|---|
| 251 | //----------------------------------------------------------------------
|
|---|
| 252 | /// Add to the vector tile_union the tiles that are in the neighbourhood
|
|---|
| 253 | /// of the specified tile_index, including itself -- start adding
|
|---|
| 254 | /// from position n_near_tiles-1, and increase n_near_tiles as
|
|---|
| 255 | /// you go along (could have done it more C++ like with vector with reserved
|
|---|
| 256 | /// space, but fear is that it would have been slower, e.g. checking
|
|---|
| 257 | /// for end of vector at each stage to decide whether to resize it)
|
|---|
| 258 | void LazyTiling9SeparateGhosts::_add_neighbours_to_tile_union(const int tile_index,
|
|---|
| 259 | vector<int> & tile_union, int & n_near_tiles) const {
|
|---|
| 260 | for (Tile3 * const * near_tile = _tiles[tile_index].begin_tiles;
|
|---|
| 261 | near_tile != _tiles[tile_index].end_tiles; near_tile++){
|
|---|
| 262 | // get the tile number
|
|---|
| 263 | tile_union[n_near_tiles] = *near_tile - & _tiles[0];
|
|---|
| 264 | n_near_tiles++;
|
|---|
| 265 | }
|
|---|
| 266 | }
|
|---|
| 267 |
|
|---|
| 268 |
|
|---|
| 269 | //----------------------------------------------------------------------
|
|---|
| 270 | /// Like _add_neighbours_to_tile_union, but only adds neighbours if
|
|---|
| 271 | /// their "tagged" status is false; when a neighbour is added its
|
|---|
| 272 | /// tagged status is set to true.
|
|---|
| 273 | inline void LazyTiling9SeparateGhosts::_add_untagged_neighbours_to_tile_union(
|
|---|
| 274 | const int tile_index,
|
|---|
| 275 | vector<int> & tile_union, int & n_near_tiles) {
|
|---|
| 276 | for (Tile3 ** near_tile = _tiles[tile_index].begin_tiles;
|
|---|
| 277 | near_tile != _tiles[tile_index].end_tiles; near_tile++){
|
|---|
| 278 | if (! (*near_tile)->tagged) {
|
|---|
| 279 | (*near_tile)->tagged = true;
|
|---|
| 280 | // get the tile number
|
|---|
| 281 | tile_union[n_near_tiles] = *near_tile - & _tiles[0];
|
|---|
| 282 | n_near_tiles++;
|
|---|
| 283 | }
|
|---|
| 284 | }
|
|---|
| 285 | }
|
|---|
| 286 |
|
|---|
| 287 | //----------------------------------------------------------------------
|
|---|
| 288 | /// Like _add_neighbours_to_tile_union, but adds tiles that are
|
|---|
| 289 | /// "neighbours" of a jet (rather than a tile) and only if a
|
|---|
| 290 | /// neighbouring tile's max_NN_dist is >= the distance between the jet
|
|---|
| 291 | /// and the nearest point on the tile. It ignores tiles that have
|
|---|
| 292 | /// already been tagged.
|
|---|
| 293 | inline void LazyTiling9SeparateGhosts::_add_untagged_neighbours_to_tile_union_using_max_info(
|
|---|
| 294 | const TiledJet3 * jet,
|
|---|
| 295 | vector<int> & tile_union, int & n_near_tiles) {
|
|---|
| 296 | Tile3 & tile = _tiles[jet->tile_index];
|
|---|
| 297 |
|
|---|
| 298 | for (Tile3 ** near_tile = tile.begin_tiles; near_tile != tile.end_tiles; near_tile++){
|
|---|
| 299 | if ((*near_tile)->tagged) continue;
|
|---|
| 300 | // here we are not allowed to miss a tile due to some rounding
|
|---|
| 301 | // error. We therefore allow for a margin of security
|
|---|
| 302 | double dist = _distance_to_tile(jet, *near_tile) - tile_edge_security_margin;
|
|---|
| 303 | // cout << " max info looked at tile " << *near_tile - &_tiles[0]
|
|---|
| 304 | // << ", dist = " << dist << " " << (*near_tile)->max_NN_dist
|
|---|
| 305 | // << endl;
|
|---|
| 306 | if (dist > (*near_tile)->max_NN_dist) continue;
|
|---|
| 307 |
|
|---|
| 308 | // cout << " max info tagged tile " << *near_tile - &_tiles[0] << endl;
|
|---|
| 309 | (*near_tile)->tagged = true;
|
|---|
| 310 | // get the tile number
|
|---|
| 311 | tile_union[n_near_tiles] = *near_tile - & _tiles[0];
|
|---|
| 312 | n_near_tiles++;
|
|---|
| 313 | }
|
|---|
| 314 | }
|
|---|
| 315 |
|
|---|
| 316 | ////--------TMPTMPTMPTMPTMP-----GPS TEMP--------------------
|
|---|
| 317 | //ostream & operator<<(ostream & ostr, const TiledJet3 & jet) {
|
|---|
| 318 | // ostr << "j" << setw(3) << jet._jets_index << ":pt2,rap,phi=" ; ostr.flush();
|
|---|
| 319 | // ostr << jet.kt2 << ","; ostr.flush();
|
|---|
| 320 | // ostr << jet.eta << ","; ostr.flush();
|
|---|
| 321 | // ostr << jet.phi; ostr.flush();
|
|---|
| 322 | // ostr << ", tile=" << jet.tile_index; ostr.flush();
|
|---|
| 323 | // return ostr;
|
|---|
| 324 | //}
|
|---|
| 325 |
|
|---|
| 326 |
|
|---|
| 327 | //----------------------------------------------------------------------
|
|---|
| 328 | /// returns a particle's distance to the edge of the specified tile
|
|---|
| 329 | inline double LazyTiling9SeparateGhosts::_distance_to_tile(const TiledJet3 * bj, const Tile3 * tile) const {
|
|---|
| 330 |
|
|---|
| 331 | // Note the careful way of checking the minimum potential deta:
|
|---|
| 332 | // unlike the phi case below, we don't calculate the distance to the
|
|---|
| 333 | // centre and subtract spacing/2. This is because of issue of
|
|---|
| 334 | // boundary tiles, which can extend far beyond spacing/2 in eta.
|
|---|
| 335 | // Using the positions of tile centers should instead be safe.
|
|---|
| 336 | double deta;
|
|---|
| 337 | if (_tiles[bj->tile_index].eta_centre == tile->eta_centre) deta = 0;
|
|---|
| 338 | //else deta = std::abs(bj->eta - tile->eta_centre) - 0.5*_tile_size_eta;
|
|---|
| 339 | else deta = std::abs(bj->eta - tile->eta_centre) - _tile_half_size_eta;
|
|---|
| 340 | // ------
|
|---|
| 341 | // |
|
|---|
| 342 | // A | B
|
|---|
| 343 | // ------
|
|---|
| 344 | // |
|
|---|
| 345 | // C | D
|
|---|
| 346 | // ------
|
|---|
| 347 |
|
|---|
| 348 | double dphi = std::abs(bj->phi - tile->phi_centre);
|
|---|
| 349 | if (dphi > pi) dphi = twopi-dphi;
|
|---|
| 350 | dphi -= _tile_half_size_phi;
|
|---|
| 351 | //dphi -= 0.5*_tile_size_phi;
|
|---|
| 352 | if (dphi < 0) dphi = 0;
|
|---|
| 353 |
|
|---|
| 354 | return dphi*dphi + deta*deta;
|
|---|
| 355 | }
|
|---|
| 356 |
|
|---|
| 357 |
|
|---|
| 358 |
|
|---|
| 359 |
|
|---|
| 360 | //----------------------------------------------------------------------
|
|---|
| 361 | /// looks at distance between jetX and jetI and updates the NN
|
|---|
| 362 | /// information if relevant; also pushes identity of jetI onto
|
|---|
| 363 | /// the vector of jets for minheap, to signal that it will have
|
|---|
| 364 | /// to be handled later.
|
|---|
| 365 | ///
|
|---|
| 366 | inline void LazyTiling9SeparateGhosts::_update_jetX_jetI_NN(TiledJet3 * jetX, TiledJet3 * jetI, vector<TiledJet3 *> & jets_for_minheap) {
|
|---|
| 367 | assert(! (jetX->is_ghost || jetI->is_ghost));
|
|---|
| 368 | double dist = _bj_dist(jetI,jetX);
|
|---|
| 369 | if (dist < jetI->NN_dist) {
|
|---|
| 370 | if (jetI != jetX) {
|
|---|
| 371 | jetI->NN_dist = dist;
|
|---|
| 372 | jetI->NN = jetX;
|
|---|
| 373 | // label jetI as needing heap action...
|
|---|
| 374 | if (!jetI->minheap_update_needed()) {
|
|---|
| 375 | jetI->label_minheap_update_needed();
|
|---|
| 376 | jets_for_minheap.push_back(jetI);
|
|---|
| 377 | }
|
|---|
| 378 | }
|
|---|
| 379 | }
|
|---|
| 380 | if (dist < jetX->NN_dist) {
|
|---|
| 381 | if (jetI != jetX) {
|
|---|
| 382 | jetX->NN_dist = dist;
|
|---|
| 383 | jetX->NN = jetI;}
|
|---|
| 384 | }
|
|---|
| 385 | }
|
|---|
| 386 |
|
|---|
| 387 |
|
|---|
| 388 | inline void LazyTiling9SeparateGhosts::_set_NN(TiledJet3 * jetI,
|
|---|
| 389 | vector<TiledJet3 *> & jets_for_minheap) {
|
|---|
| 390 | assert(! jetI->is_ghost);
|
|---|
| 391 | jetI->NN_dist = _R2;
|
|---|
| 392 | jetI->NN = NULL;
|
|---|
| 393 | // label jetI as needing heap action...
|
|---|
| 394 | if (!jetI->minheap_update_needed()) {
|
|---|
| 395 | jetI->label_minheap_update_needed();
|
|---|
| 396 | jets_for_minheap.push_back(jetI);}
|
|---|
| 397 | // now go over tiles that are neighbours of I (include own tile)
|
|---|
| 398 | Tile3 * tile_ptr = &_tiles[jetI->tile_index];
|
|---|
| 399 | //if (tile_ptr->is_near_zero_phi(_tile_size_phi)) {
|
|---|
| 400 | for (Tile3 ** near_tile = tile_ptr->begin_tiles;
|
|---|
| 401 | near_tile != tile_ptr->end_tiles; near_tile++) {
|
|---|
| 402 | // for own tile, this will be zero automatically: should we be clever
|
|---|
| 403 | // and skip the test? (With some doubling of code?)
|
|---|
| 404 | if (jetI->NN_dist < _distance_to_tile(jetI, *near_tile)) continue;
|
|---|
| 405 | // and then over the contents of that tile
|
|---|
| 406 | for (TiledJet3 * jetJ = (*near_tile)->head; jetJ != NULL; jetJ = jetJ->next) {
|
|---|
| 407 | double dist = _bj_dist(jetI,jetJ);
|
|---|
| 408 | if (dist < jetI->NN_dist && jetJ != jetI) {
|
|---|
| 409 | jetI->NN_dist = dist; jetI->NN = jetJ;
|
|---|
| 410 | }
|
|---|
| 411 | }
|
|---|
| 412 | // deal with the ghosts
|
|---|
| 413 | for (TiledJet3 * jetJ = (*near_tile)->ghost_head; jetJ != NULL; jetJ = jetJ->next) {
|
|---|
| 414 | double dist = _bj_dist(jetI,jetJ);
|
|---|
| 415 | if (dist < jetI->NN_dist) {
|
|---|
| 416 | jetI->NN_dist = dist; jetI->NN = jetJ;
|
|---|
| 417 | }
|
|---|
| 418 | }
|
|---|
| 419 | }
|
|---|
| 420 | // } else {
|
|---|
| 421 | // // second copy that exploits the fact that for this tile we needn't worry
|
|---|
| 422 | // // about periodicity
|
|---|
| 423 | // for (Tile3 ** near_tile = tile_ptr->begin_tiles;
|
|---|
| 424 | // near_tile != tile_ptr->end_tiles; near_tile++) {
|
|---|
| 425 | // // for own tile, this will be zero automatically: should we be clever
|
|---|
| 426 | // // and skip the test? (With some doubling of code?)
|
|---|
| 427 | // if (jetI->NN_dist < _distance_to_tile(jetI, *near_tile)) continue;
|
|---|
| 428 | // // and then over the contents of that tile
|
|---|
| 429 | // for (TiledJet3 * jetJ = (*near_tile)->head;
|
|---|
| 430 | // jetJ != NULL; jetJ = jetJ->next) {
|
|---|
| 431 | // double dist = _bj_dist_not_periodic(jetI,jetJ);
|
|---|
| 432 | // if (dist < jetI->NN_dist && jetJ != jetI) {
|
|---|
| 433 | // jetI->NN_dist = dist; jetI->NN = jetJ;
|
|---|
| 434 | // }
|
|---|
| 435 | // }
|
|---|
| 436 | // }
|
|---|
| 437 | // }
|
|---|
| 438 | }
|
|---|
| 439 |
|
|---|
| 440 |
|
|---|
| 441 | void LazyTiling9SeparateGhosts::run() {
|
|---|
| 442 |
|
|---|
| 443 | //_initialise_tiles();
|
|---|
| 444 |
|
|---|
| 445 | int ntot = _jets.size();
|
|---|
| 446 | if (ntot == 0) return;
|
|---|
| 447 |
|
|---|
| 448 | TiledJet3 * briefjets = new TiledJet3[ntot];
|
|---|
| 449 | TiledJet3 * jetA = briefjets, * jetB;
|
|---|
| 450 | // avoid warning about uninitialised oldB below;
|
|---|
| 451 | // only valid for ntot>=1 (hence the test ntot==0 test above)
|
|---|
| 452 | TiledJet3 oldB = briefjets[0];
|
|---|
| 453 |
|
|---|
| 454 |
|
|---|
| 455 | // will be used quite deep inside loops, but declare it here so that
|
|---|
| 456 | // memory (de)allocation gets done only once
|
|---|
| 457 | vector<int> tile_union(3*n_tile_neighbours);
|
|---|
| 458 |
|
|---|
| 459 | TiledJet3 * head = briefjets; // a nicer way of naming start
|
|---|
| 460 |
|
|---|
| 461 | // initialise the basic jet info
|
|---|
| 462 | //
|
|---|
| 463 | // Note that the threshold is a static member of the class
|
|---|
| 464 | // first get the particles we'll keep as "real"
|
|---|
| 465 | for (int i = 0; i< ntot; i++) {
|
|---|
| 466 | bool is_ghost = _jets[i].perp2() < ghost_pt2_threshold;
|
|---|
| 467 | if (!is_ghost) {
|
|---|
| 468 | _tj_set_jetinfo(jetA, i, is_ghost);
|
|---|
| 469 | jetA++; // move on to next entry of briefjets
|
|---|
| 470 | }
|
|---|
| 471 | }
|
|---|
| 472 | int nreal = jetA - briefjets;
|
|---|
| 473 | // then the ones we will label as ghosts
|
|---|
| 474 | for (int i = 0; i< ntot; i++) {
|
|---|
| 475 | bool is_ghost = _jets[i].perp2() < ghost_pt2_threshold;
|
|---|
| 476 | if (is_ghost) {
|
|---|
| 477 | _tj_set_jetinfo(jetA, i, is_ghost);
|
|---|
| 478 | jetA++; // move on to next entry of briefjets
|
|---|
| 479 | }
|
|---|
| 480 | }
|
|---|
| 481 |
|
|---|
| 482 |
|
|---|
| 483 | // set up the initial nearest neighbour information
|
|---|
| 484 | vector<Tile3>::iterator tile;
|
|---|
| 485 | for (tile = _tiles.begin(); tile != _tiles.end(); tile++) {
|
|---|
| 486 | // first do it on this tile
|
|---|
| 487 | for (jetA = tile->head; jetA != NULL; jetA = jetA->next) {
|
|---|
| 488 | // real particles cluster with real particles
|
|---|
| 489 | for (jetB = tile->head; jetB != jetA; jetB = jetB->next) {
|
|---|
| 490 | double dist = _bj_dist_not_periodic(jetA,jetB);
|
|---|
| 491 | if (dist < jetA->NN_dist) {jetA->NN_dist = dist; jetA->NN = jetB;}
|
|---|
| 492 | if (dist < jetB->NN_dist) {jetB->NN_dist = dist; jetB->NN = jetA;}
|
|---|
| 493 | }
|
|---|
| 494 | // they can also cluster with ghosts
|
|---|
| 495 | for (jetB = tile->ghost_head; jetB != NULL; jetB = jetB->next) {
|
|---|
| 496 | double dist = _bj_dist_not_periodic(jetA,jetB);
|
|---|
| 497 | if (dist < jetA->NN_dist) {jetA->NN_dist = dist; jetA->NN = jetB;}
|
|---|
| 498 | }
|
|---|
| 499 | }
|
|---|
| 500 | // only look out for NN dists of real particles, because ghosts never
|
|---|
| 501 | // have NN in our structure
|
|---|
| 502 | for (jetA = tile->head; jetA != NULL; jetA = jetA->next) {
|
|---|
| 503 | if (jetA->NN_dist > tile->max_NN_dist) tile->max_NN_dist = jetA->NN_dist;
|
|---|
| 504 | }
|
|---|
| 505 | }
|
|---|
| 506 | for (tile = _tiles.begin(); tile != _tiles.end(); tile++) {
|
|---|
| 507 | //if (tile->is_near_zero_phi(_tile_size_phi)) {
|
|---|
| 508 | // then do it for RH tiles;
|
|---|
| 509 | for (Tile3 ** RTile = tile->RH_tiles; RTile != tile->end_tiles; RTile++) {
|
|---|
| 510 | for (jetA = tile->head; jetA != NULL; jetA = jetA->next) {
|
|---|
| 511 | double dist_to_tile = _distance_to_tile(jetA, *RTile);
|
|---|
| 512 | // it only makes sense to do a tile if jetA is close enough to the Rtile
|
|---|
| 513 | // either for a jet in the Rtile to be closer to jetA than it's current NN
|
|---|
| 514 | // or if jetA could be closer to something in the Rtile than the largest
|
|---|
| 515 | // NN distance within the RTile.
|
|---|
| 516 | //
|
|---|
| 517 | // GPS note: also tried approach where we perform only the
|
|---|
| 518 | // first test and run over all surrounding tiles
|
|---|
| 519 | // (not just RH ones). The test is passed less
|
|---|
| 520 | // frequently, but one is running over more tiles
|
|---|
| 521 | // and on balance, for the trial event we used, it's
|
|---|
| 522 | // a bit slower.
|
|---|
| 523 | bool relevant_for_jetA = dist_to_tile <= jetA->NN_dist;
|
|---|
| 524 | bool relevant_for_RTile = dist_to_tile <= (*RTile)->max_NN_dist;
|
|---|
| 525 | if (relevant_for_jetA || relevant_for_RTile) {
|
|---|
| 526 | for (jetB = (*RTile)->head; jetB != NULL; jetB = jetB->next) {
|
|---|
| 527 | double dist = _bj_dist(jetA,jetB);
|
|---|
| 528 | if (dist < jetA->NN_dist) {jetA->NN_dist = dist; jetA->NN = jetB;}
|
|---|
| 529 | if (dist < jetB->NN_dist) {jetB->NN_dist = dist; jetB->NN = jetA;}
|
|---|
| 530 | }
|
|---|
| 531 | }
|
|---|
| 532 | // now do the check over ghosts
|
|---|
| 533 | if (relevant_for_jetA) {
|
|---|
| 534 | for (jetB = (*RTile)->ghost_head; jetB != NULL; jetB = jetB->next) {
|
|---|
| 535 | double dist = _bj_dist(jetA,jetB);
|
|---|
| 536 | if (dist < jetA->NN_dist) {jetA->NN_dist = dist; jetA->NN = jetB;}
|
|---|
| 537 | }
|
|---|
| 538 | }
|
|---|
| 539 | }
|
|---|
| 540 | }
|
|---|
| 541 | // and do a special loop to catch ghosts that are among the LH tiles
|
|---|
| 542 | for (Tile3 ** LTile = tile->surrounding_tiles; LTile != tile->RH_tiles; LTile++) {
|
|---|
| 543 | for (jetA = tile->head; jetA != NULL; jetA = jetA->next) {
|
|---|
| 544 | double dist_to_tile = _distance_to_tile(jetA, *LTile);
|
|---|
| 545 | // it only makes sense to do a tile if jetA is close enough to the Rtile
|
|---|
| 546 | // for a (ghost) jet in the Ltile to be closer to jetA than it's current NN.
|
|---|
| 547 | bool relevant_for_jetA = dist_to_tile <= jetA->NN_dist;
|
|---|
| 548 | if (relevant_for_jetA) {
|
|---|
| 549 | for (jetB = (*LTile)->ghost_head; jetB != NULL; jetB = jetB->next) {
|
|---|
| 550 | double dist = _bj_dist(jetA,jetB);
|
|---|
| 551 | if (dist < jetA->NN_dist) {jetA->NN_dist = dist; jetA->NN = jetB;}
|
|---|
| 552 | }
|
|---|
| 553 | }
|
|---|
| 554 | }
|
|---|
| 555 | }
|
|---|
| 556 | // } else {
|
|---|
| 557 | // // this second version of the code uses the faster
|
|---|
| 558 | // // "not_periodic" version because it knows that the tile is
|
|---|
| 559 | // // sufficiently far from the edge.
|
|---|
| 560 | // for (Tile3 ** RTile = tile->RH_tiles; RTile != tile->end_tiles; RTile++) {
|
|---|
| 561 | // for (jetA = tile->head; jetA != NULL; jetA = jetA->next) {
|
|---|
| 562 | // double dist_to_tile = _distance_to_tile(jetA, *RTile);
|
|---|
| 563 | // bool relevant_for_jetA = dist_to_tile <= jetA->NN_dist;
|
|---|
| 564 | // bool relevant_for_RTile = dist_to_tile <= (*RTile)->max_NN_dist;
|
|---|
| 565 | // if (relevant_for_jetA || relevant_for_RTile) {
|
|---|
| 566 | // for (jetB = (*RTile)->head; jetB != NULL; jetB = jetB->next) {
|
|---|
| 567 | // double dist = _bj_dist_not_periodic(jetA,jetB);
|
|---|
| 568 | // if (dist < jetA->NN_dist) {jetA->NN_dist = dist; jetA->NN = jetB;}
|
|---|
| 569 | // if (dist < jetB->NN_dist) {jetB->NN_dist = dist; jetB->NN = jetA;}
|
|---|
| 570 | // }
|
|---|
| 571 | // }
|
|---|
| 572 | // }
|
|---|
| 573 | // }
|
|---|
| 574 | // }
|
|---|
| 575 | // no need to do it for LH tiles, since they are implicitly done
|
|---|
| 576 | // when we set NN for both jetA and jetB on the RH tiles.
|
|---|
| 577 | }
|
|---|
| 578 | // Now update the max_NN_dist within each tile. Not strictly
|
|---|
| 579 | // necessary, because existing max_NN_dist is an upper bound. but
|
|---|
| 580 | // costs little and may give some efficiency gain later.
|
|---|
| 581 | // (Do it only for real particles -- ghosts don't come into the game).
|
|---|
| 582 | for (tile = _tiles.begin(); tile != _tiles.end(); tile++) {
|
|---|
| 583 | tile->max_NN_dist = 0;
|
|---|
| 584 | for (jetA = tile->head; jetA != NULL; jetA = jetA->next) {
|
|---|
| 585 | if (jetA->NN_dist > tile->max_NN_dist) tile->max_NN_dist = jetA->NN_dist;
|
|---|
| 586 | }
|
|---|
| 587 | }
|
|---|
| 588 |
|
|---|
| 589 |
|
|---|
| 590 | vector<double> diJs(nreal);
|
|---|
| 591 | for (int i = 0; i < nreal; i++) {
|
|---|
| 592 | diJs[i] = _bj_diJ(&briefjets[i]);
|
|---|
| 593 | briefjets[i].label_minheap_update_done();
|
|---|
| 594 | }
|
|---|
| 595 | MinHeap minheap(diJs);
|
|---|
| 596 | // have a stack telling us which jets we'll have to update on the heap
|
|---|
| 597 | vector<TiledJet3 *> jets_for_minheap;
|
|---|
| 598 | jets_for_minheap.reserve(ntot); // GPS TMP: return here
|
|---|
| 599 |
|
|---|
| 600 | // now run the recombination loop
|
|---|
| 601 | while (nreal > 0) {
|
|---|
| 602 |
|
|---|
| 603 | double diJ_min = minheap.minval() *_invR2;
|
|---|
| 604 | jetA = head + minheap.minloc();
|
|---|
| 605 |
|
|---|
| 606 | // do the recombination between A and B
|
|---|
| 607 | jetB = jetA->NN;
|
|---|
| 608 |
|
|---|
| 609 | if (jetB != NULL) {
|
|---|
| 610 | // jet-jet recombination
|
|---|
| 611 | // If necessary relabel A & B to ensure jetB < jetA, that way if
|
|---|
| 612 | // the larger of them == newtail then that ends up being jetA and
|
|---|
| 613 | // the new jet that is added as jetB is inserted in a position that
|
|---|
| 614 | // has a future!
|
|---|
| 615 | if (jetA < jetB) {std::swap(jetA,jetB);}
|
|---|
| 616 |
|
|---|
| 617 | int nn; // new jet index
|
|---|
| 618 | _cs.plugin_record_ij_recombination(jetA->_jets_index, jetB->_jets_index, diJ_min, nn);
|
|---|
| 619 |
|
|---|
| 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 | bool is_ghost = false;
|
|---|
| 625 | _tj_set_jetinfo(jetB, nn, is_ghost); // cause jetB to become _jets[nn]
|
|---|
| 626 | // (also registers the jet in the tiling); it will never be a ghost
|
|---|
| 627 | } else {
|
|---|
| 628 | // jet-beam recombination
|
|---|
| 629 | // get the hist_index
|
|---|
| 630 | _cs.plugin_record_iB_recombination(jetA->_jets_index, diJ_min);
|
|---|
| 631 | _bj_remove_from_tiles(jetA);
|
|---|
| 632 | }
|
|---|
| 633 |
|
|---|
| 634 | // remove the minheap entry for jetA if jetA is a real particle;
|
|---|
| 635 | // only in this case
|
|---|
| 636 | if (!jetA->is_ghost) {
|
|---|
| 637 | minheap.remove(jetA-head);
|
|---|
| 638 | // jetB cannot be a ghost; so when jetA is not a ghost we decrease
|
|---|
| 639 | // our count of remaining real particles
|
|---|
| 640 | nreal--;
|
|---|
| 641 | }
|
|---|
| 642 |
|
|---|
| 643 |
|
|---|
| 644 |
|
|---|
| 645 |
|
|---|
| 646 | // first establish the set of tiles over which we are going to
|
|---|
| 647 | // have to run searches for updated and new nearest-neighbours --
|
|---|
| 648 | // basically a combination of vicinity of the tiles of the two old
|
|---|
| 649 | // and one new jet.
|
|---|
| 650 | int n_near_tiles = 0;
|
|---|
| 651 | // add tiles neighbouring A even if A is a ghost because a ghost
|
|---|
| 652 | // could have been the NN of some other particle
|
|---|
| 653 | _add_untagged_neighbours_to_tile_union_using_max_info(jetA,
|
|---|
| 654 | tile_union, n_near_tiles);
|
|---|
| 655 | if (jetB != NULL) {
|
|---|
| 656 | _add_untagged_neighbours_to_tile_union_using_max_info(&oldB,
|
|---|
| 657 | tile_union,n_near_tiles);
|
|---|
| 658 | jetB->label_minheap_update_needed();
|
|---|
| 659 | jets_for_minheap.push_back(jetB);
|
|---|
| 660 | }
|
|---|
| 661 |
|
|---|
| 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 | if (jetB != NULL) {
|
|---|
| 667 | Tile3 & jetB_tile = _tiles[jetB->tile_index];
|
|---|
| 668 | for (Tile3 ** near_tile = jetB_tile.begin_tiles;
|
|---|
| 669 | near_tile != jetB_tile.end_tiles; near_tile++) {
|
|---|
| 670 |
|
|---|
| 671 | double dist_to_tile = _distance_to_tile(jetB, *near_tile);
|
|---|
| 672 | // use <= in next line so that on first tile, relevant_for_jetB is
|
|---|
| 673 | // set to true
|
|---|
| 674 | bool relevant_for_jetB = dist_to_tile <= jetB->NN_dist;
|
|---|
| 675 | bool relevant_for_near_tile = dist_to_tile <= (*near_tile)->max_NN_dist;
|
|---|
| 676 | bool relevant = relevant_for_jetB || relevant_for_near_tile;
|
|---|
| 677 | // this first option decides exactly what loop to do based on whether
|
|---|
| 678 | // the near tile was tagged. You'd think it's more efficient, but
|
|---|
| 679 | // not necessarily...
|
|---|
| 680 | if (relevant) {
|
|---|
| 681 | if ((*near_tile)->tagged) {
|
|---|
| 682 | for (TiledJet3 * jetI = (*near_tile)->head; jetI != NULL; jetI = jetI->next) {
|
|---|
| 683 | if (jetI->NN == jetA || jetI->NN == jetB) _set_NN(jetI, jets_for_minheap);
|
|---|
| 684 | _update_jetX_jetI_NN(jetB, jetI, jets_for_minheap);
|
|---|
| 685 | }
|
|---|
| 686 | (*near_tile)->tagged = false;
|
|---|
| 687 | } else {
|
|---|
| 688 | for (TiledJet3 * jetI = (*near_tile)->head; jetI != NULL; jetI = jetI->next) {
|
|---|
| 689 | _update_jetX_jetI_NN(jetB, jetI, jets_for_minheap);
|
|---|
| 690 | }
|
|---|
| 691 | }
|
|---|
| 692 | }
|
|---|
| 693 | // Now HANDLE GHOSTS
|
|---|
| 694 | if (relevant_for_jetB) {
|
|---|
| 695 | for (TiledJet3 * jetI = (*near_tile)->ghost_head; jetI != NULL; jetI = jetI->next) {
|
|---|
| 696 | double dist = _bj_dist(jetB,jetI);
|
|---|
| 697 | if (dist < jetB->NN_dist) {
|
|---|
| 698 | jetB->NN_dist = dist; jetB->NN = jetI;
|
|---|
| 699 | }
|
|---|
| 700 | }
|
|---|
| 701 | }
|
|---|
| 702 |
|
|---|
| 703 |
|
|---|
| 704 | // this second option does everything independently of whether the near tile
|
|---|
| 705 | // was tagged -- somehow you'd expect it to be slower, but it may actually be
|
|---|
| 706 | // marginally faster.
|
|---|
| 707 | // if (relevant_for_jetB || relevant_for_near_tile) {
|
|---|
| 708 | // for (TiledJet3 * jetI = (*near_tile)->head; jetI != NULL; jetI = jetI->next) {
|
|---|
| 709 | //
|
|---|
| 710 | // if (jetI->NN == jetA || (jetI->NN == jetB && jetB != NULL)) {
|
|---|
| 711 | // _set_NN(jetI, jets_for_minheap);
|
|---|
| 712 | // }
|
|---|
| 713 | //
|
|---|
| 714 | // _update_jetX_jetI_NN(jetB, jetI, jets_for_minheap);
|
|---|
| 715 | // // -- Keep this old inline code for later speed tests
|
|---|
| 716 | // // double dist = _bj_dist(jetI,jetB);
|
|---|
| 717 | // // if (dist < jetI->NN_dist) {
|
|---|
| 718 | // // if (jetI != jetB) {
|
|---|
| 719 | // // jetI->NN_dist = dist;
|
|---|
| 720 | // // jetI->NN = jetB;
|
|---|
| 721 | // // // label jetI as needing heap action...
|
|---|
| 722 | // // if (!jetI->minheap_update_needed()) {
|
|---|
| 723 | // // jetI->label_minheap_update_needed();
|
|---|
| 724 | // // jets_for_minheap.push_back(jetI);
|
|---|
| 725 | // // }
|
|---|
| 726 | // // }
|
|---|
| 727 | // // }
|
|---|
| 728 | // // if (dist < jetB->NN_dist) {
|
|---|
| 729 | // // if (jetI != jetB) {
|
|---|
| 730 | // // jetB->NN_dist = dist;
|
|---|
| 731 | // // jetB->NN = jetI;}
|
|---|
| 732 | // // }
|
|---|
| 733 | // }
|
|---|
| 734 | // (*near_tile)->tagged = false;
|
|---|
| 735 | // }
|
|---|
| 736 | }
|
|---|
| 737 | }
|
|---|
| 738 |
|
|---|
| 739 | // now run over the tiles that were tagged earlier and that we haven't yet
|
|---|
| 740 | // had a change to visit.
|
|---|
| 741 | for (int itile = 0; itile < n_near_tiles; itile++) {
|
|---|
| 742 | Tile3 * tile_ptr = &_tiles[tile_union[itile]];
|
|---|
| 743 | if (!tile_ptr->tagged) continue; // because earlier loop may have undone the tag
|
|---|
| 744 | tile_ptr->tagged = false;
|
|---|
| 745 | // run over all jets in the current tile
|
|---|
| 746 | for (TiledJet3 * jetI = tile_ptr->head; jetI != NULL; jetI = jetI->next) {
|
|---|
| 747 | // see if jetI had jetA or jetB as a NN -- if so recalculate the NN
|
|---|
| 748 | if (jetI->NN == jetA || (jetI->NN == jetB && jetB != NULL)) {
|
|---|
| 749 | _set_NN(jetI, jets_for_minheap);
|
|---|
| 750 | }
|
|---|
| 751 | }
|
|---|
| 752 | }
|
|---|
| 753 |
|
|---|
| 754 | // deal with jets whose minheap entry needs updating
|
|---|
| 755 | //if (verbose) cout << " jets whose NN was modified: " << endl;
|
|---|
| 756 | while (jets_for_minheap.size() > 0) {
|
|---|
| 757 | TiledJet3 * jetI = jets_for_minheap.back();
|
|---|
| 758 | jets_for_minheap.pop_back();
|
|---|
| 759 | minheap.update(jetI-head, _bj_diJ(jetI));
|
|---|
| 760 | jetI->label_minheap_update_done();
|
|---|
| 761 | // handle max_NN_dist update for all jets that might have
|
|---|
| 762 | // seen a change (increase) of distance
|
|---|
| 763 | Tile3 & tile_I = _tiles[jetI->tile_index];
|
|---|
| 764 | if (tile_I.max_NN_dist < jetI->NN_dist) tile_I.max_NN_dist = jetI->NN_dist;
|
|---|
| 765 | }
|
|---|
| 766 | }
|
|---|
| 767 |
|
|---|
| 768 | // final cleaning up;
|
|---|
| 769 | delete[] briefjets;
|
|---|
| 770 | }
|
|---|
| 771 |
|
|---|
| 772 |
|
|---|
| 773 | FASTJET_END_NAMESPACE
|
|---|
