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