LazyTiling25.hh

#ifndef __FASTJET_LAZYTILING25_HH__
#define __FASTJET_LAZYTILING25_HH__

// #define INSTRUMENT2 1

//FJSTARTHEADER
// $Id: LazyTiling25.hh 3477 2014-07-29 14:34:39Z salam $
//
// Copyright (c) 2005-2014, 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"
#include "fastjet/internal/LazyTiling9Alt.hh"
#include "fastjet/internal/LazyTiling9.hh"



FASTJET_BEGIN_NAMESPACE      // defined in fastjet/internal/base.hh

typedef Tile2Base<25> Tile25;


// class Tile25 {
// public:
//   /// pointers to neighbouring tiles, including self
//   Tile25 *   begin_tiles[25];
//   /// neighbouring tiles, excluding self
//   Tile25 **  surrounding_tiles; 
//   /// half of neighbouring tiles, no self
//   Tile25 **  RH_tiles;  
//   /// just beyond end of tiles
//   Tile25 **  end_tiles; 
//   /// start of list of BriefJets contained in this tile
//   TiledJet * head;    
//   /// sometimes useful to be able to tag a tile
//   bool     tagged;    
//   /// for all particles in the tile, this stores the largest of the
//   /// (squared) nearest-neighbour distances.
//   double max_NN_dist;
//   double eta_centre, phi_centre;
// 
//   bool is_near_zero_phi(double tile_size_phi) const {
//     return phi_centre < 2*tile_size_phi || (twopi-phi_centre) < 2*tile_size_phi;
//   }
// };


//----------------------------------------------------------------------
class LazyTiling25 {
public:
  LazyTiling25(ClusterSequence & cs);

  void run();

protected:
  ClusterSequence & _cs;
  const std::vector<PseudoJet> & _jets;
  std::vector<Tile25> _tiles;

#ifdef INSTRUMENT2
  int _ncall; // GPS tmp
  int _ncall_dtt; // GPS tmp
#endif // INSTRUMENT2

  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 Tile25 *) 
#ifdef INSTRUMENT2
    ;
#else
    const;
#endif 
  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) 
#ifdef INSTRUMENT2
    {
    _ncall++; // GPS tmp
#else
    const {
#endif 
    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)
#ifdef INSTRUMENT2
    {
    _ncall++; // GPS tmp
#else
    const {
#endif 
    //_ncall++; // GPS tmp
    double dphi = jetA->phi - jetB->phi;
    double deta = (jetA->eta - jetB->eta);
    return dphi*dphi + deta*deta;
  }

};


FASTJET_END_NAMESPACE

#endif // __FASTJET_LAZYTILING25_HH__