Loading...
Searching...
No Matches
LazyTiling9Alt.hh
#ifndef __FASTJET_LAZYTILING9ALT_HH__
#define __FASTJET_LAZYTILING9ALT_HH__
//FJSTARTHEADER
// $Id$
//
// Copyright (c) 2005-2021, Matteo Cacciari, Gavin P. Salam and Gregory Soyez
//
//----------------------------------------------------------------------
// This file is part of FastJet.
//
// FastJet is free software; you can redistribute it and/or modify
// it under the terms of the GNU General Public License as published by
// the Free Software Foundation; either version 2 of the License, or
// (at your option) any later version.
//
// The algorithms that underlie FastJet have required considerable
// development. They are described in the original FastJet paper,
// hep-ph/0512210 and in the manual, arXiv:1111.6097. If you use
// FastJet as part of work towards a scientific publication, please
// quote the version you use and include a citation to the manual and
// optionally also to hep-ph/0512210.
//
// FastJet is distributed in the hope that it will be useful,
// but WITHOUT ANY WARRANTY; without even the implied warranty of
// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
// GNU General Public License for more details.
//
// You should have received a copy of the GNU General Public License
// along with FastJet. If not, see <http://www.gnu.org/licenses/>.
//----------------------------------------------------------------------
//FJENDHEADER
//#include "fastjet/PseudoJet.hh"
#include "fastjet/internal/MinHeap.hh"
#include "fastjet/ClusterSequence.hh"
FASTJET_BEGIN_NAMESPACE // defined in fastjet/internal/base.hh
/// Rounding errors in the Lazy strategies may cause the following
/// problem: when browsing tiles in the vicinity of the particles
/// being clustered in order to decide which of these tiles may
/// contain particles that need to be updated (because theit NN is one
/// of the particles that are currently clustered), we discard tiles
/// that are deemed "too far from the cell" by the "max_NN_dist"
/// criterion. Because of rounding error, this condition can sometimes
/// miss cases where an update is needed.
///
/// An example of this happens if a particle '1' is, say, at the lower
/// edge of the rapidity of a given tile, with a particle '2' in the
/// tile directly on its left at the same rapidity. Assume also that
/// max_NN_dist in 2's tile corresponds to the distance between 2 and
/// teh tile of 1. If 2 is 1's NN then in case 2 gets clustered, 1's
/// NN needs to be updated. However, rounding errors in the
/// calculation of the distance between 1 and 2 may result is
/// something slightly larger than the max_NN_dist in 2's tile.
///
/// This situation corresponds to the bug reported by Jochen Olt on
/// February 12 2015 [see issue-tracker/2015-02-infinite-loop],
/// causing an infinite loop.
///
/// To prevent this, the simplest solution is, when looking at tiles
/// to browse for updateds, to add a margin of security close to the
/// edges of the cell, i.e. instead of updating only tiles for which
/// distance<=max_NN_dist, we will update tiles for which
/// distance<=max_NN_dist+tile_edge_security_margin.
///
/// Note that this does not need to be done when computing nearest
/// neighbours [rounding errors are tolerated there] but it is
/// critical when tracking points that have to be updated.
const double tile_edge_security_margin=1.0e-7;
/// structure analogous to BriefJet, but with the extra information
/// needed for dealing with tiles
class TiledJet {
public:
double eta, phi, kt2, NN_dist;
TiledJet * NN, *previous, * next;
int _jets_index, tile_index;
bool _minheap_update_needed;
// indicate whether jets need to have their minheap entries
// updated).
inline void label_minheap_update_needed() {_minheap_update_needed = true;}
inline void label_minheap_update_done() {_minheap_update_needed = false;}
inline bool minheap_update_needed() const {return _minheap_update_needed;}
};
const int n_tile_neighbours = 9;
class Tile {
public:
typedef double (Tile::*DistToTileFn)(const TiledJet*) const;
typedef std::pair<Tile *, DistToTileFn> TileFnPair;
/// pointers to neighbouring tiles, including self
TileFnPair begin_tiles[n_tile_neighbours];
/// neighbouring tiles, excluding self
TileFnPair * surrounding_tiles;
/// half of neighbouring tiles, no self
TileFnPair * RH_tiles;
/// just beyond end of tiles
TileFnPair * end_tiles;
/// start of list of BriefJets contained in this tile
TiledJet * head;
/// sometimes useful to be able to tag a tile
bool tagged;
/// true for tiles where the delta phi calculation needs
/// potentially to account for periodicity in phi
bool use_periodic_delta_phi;
/// for all particles in the tile, this stores the largest of the
/// (squared) nearest-neighbour distances.
double max_NN_dist;
double eta_min, eta_max, phi_min, phi_max;
double distance_to_centre(const TiledJet *) const {return 0;}
double distance_to_left(const TiledJet * jet) const {
double deta = jet->eta - eta_min;
return deta*deta;
}
double distance_to_right(const TiledJet * jet) const {
double deta = jet->eta - eta_max;
return deta*deta;
}
double distance_to_bottom(const TiledJet * jet) const {
double dphi = jet->phi - phi_min;
return dphi*dphi;
}
double distance_to_top(const TiledJet * jet) const {
double dphi = jet->phi - phi_max;
return dphi*dphi;
}
double distance_to_left_top(const TiledJet * jet) const {
double deta = jet->eta - eta_min;
double dphi = jet->phi - phi_max;
return deta*deta + dphi*dphi;
}
double distance_to_left_bottom(const TiledJet * jet) const {
double deta = jet->eta - eta_min;
double dphi = jet->phi - phi_min;
return deta*deta + dphi*dphi;
}
double distance_to_right_top(const TiledJet * jet) const {
double deta = jet->eta - eta_max;
double dphi = jet->phi - phi_max;
return deta*deta + dphi*dphi;
}
double distance_to_right_bottom(const TiledJet * jet) const {
double deta = jet->eta - eta_max;
double dphi = jet->phi - phi_min;
return deta*deta + dphi*dphi;
}
};
//----------------------------------------------------------------------
class LazyTiling9Alt {
public:
LazyTiling9Alt(ClusterSequence & cs);
void run();
//void get_next_clustering(int & jetA_index, int & jetB_index, double & dij);
protected:
ClusterSequence & _cs;
const std::vector<PseudoJet> & _jets;
std::vector<Tile> _tiles;
double _Rparam, _R2, _invR2;
double _tiles_eta_min, _tiles_eta_max;
double _tile_size_eta, _tile_size_phi;
double _tile_half_size_eta, _tile_half_size_phi;
int _n_tiles_phi,_tiles_ieta_min,_tiles_ieta_max;
std::vector<TiledJet *> _jets_for_minheap;
//MinHeap _minheap;
void _initialise_tiles();
// reasonably robust return of tile index given ieta and iphi, in particular
// it works even if iphi is negative
inline int _tile_index (int ieta, int iphi) const {
// note that (-1)%n = -1 so that we have to add _n_tiles_phi
// before performing modulo operation
return (ieta-_tiles_ieta_min)*_n_tiles_phi
+ (iphi+_n_tiles_phi) % _n_tiles_phi;
}
void _bj_remove_from_tiles(TiledJet * const jet);
/// returns the tile index given the eta and phi values of a jet
int _tile_index(const double eta, const double phi) const;
// sets up information regarding the tiling of the given jet
void _tj_set_jetinfo(TiledJet * const jet, const int _jets_index);
void _print_tiles(TiledJet * briefjets ) const;
void _add_neighbours_to_tile_union(const int tile_index,
std::vector<int> & tile_union, int & n_near_tiles) const;
void _add_untagged_neighbours_to_tile_union(const int tile_index,
std::vector<int> & tile_union, int & n_near_tiles);
void _add_untagged_neighbours_to_tile_union_using_max_info(const TiledJet * const jet,
std::vector<int> & tile_union, int & n_near_tiles);
//double _distance_to_tile(const TiledJet * bj, const Tile *) const;
void _update_jetX_jetI_NN(TiledJet * jetX, TiledJet * jetI, std::vector<TiledJet *> & jets_for_minheap);
void _set_NN(TiledJet * jetI, std::vector<TiledJet *> & jets_for_minheap);
// return the diJ (multiplied by _R2) for this jet assuming its NN
// info is correct
template <class J> double _bj_diJ(const J * const jet) const {
double kt2 = jet->kt2;
if (jet->NN != NULL) {if (jet->NN->kt2 < kt2) {kt2 = jet->NN->kt2;}}
return jet->NN_dist * kt2;
}
//----------------------------------------------------------------------
template <class J> inline void _bj_set_jetinfo(
J * const jetA, const int _jets_index) const {
jetA->eta = _jets[_jets_index].rap();
jetA->phi = _jets[_jets_index].phi_02pi();
jetA->kt2 = _cs.jet_scale_for_algorithm(_jets[_jets_index]);
jetA->_jets_index = _jets_index;
// initialise NN info as well
jetA->NN_dist = _R2;
jetA->NN = NULL;
}
//----------------------------------------------------------------------
template <class J> inline double _bj_dist(
const J * const jetA, const J * const jetB) const {
double dphi = std::abs(jetA->phi - jetB->phi);
double deta = (jetA->eta - jetB->eta);
if (dphi > pi) {dphi = twopi - dphi;}
return dphi*dphi + deta*deta;
}
//----------------------------------------------------------------------
template <class J> inline double _bj_dist_not_periodic(
const J * const jetA, const J * const jetB) const {
double dphi = jetA->phi - jetB->phi;
double deta = (jetA->eta - jetB->eta);
return dphi*dphi + deta*deta;
}
};
FASTJET_END_NAMESPACE
#endif // __FASTJET_LAZYTILING9ALT_HH__
fastjet::TiledJet
structure analogous to BriefJet, but with the extra information needed for dealing with tiles
Definition LazyTiling9Alt.hh:75
fastjet::tile_edge_security_margin
const double tile_edge_security_margin
Rounding errors in the Lazy strategies may cause the following problem: when browsing tiles in the vi...
Definition LazyTiling9Alt.hh:71