ClusterSequence_TiledN2.cc

//FJSTARTHEADER
// $Id: ClusterSequence_TiledN2.cc 4420 2019-11-29 09:28:20Z soyez $
//
// Copyright (c) 2005-2019, 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


// The tiled N^2 part of the ClusterSequence class -- separated out
// from the rest of the class implementation so as to speed up
// compilation of this particular part while it is under test.

#include<iostream>
#include<vector>
#include<cmath>
#include<algorithm>
#include "fastjet/PseudoJet.hh"
#include "fastjet/ClusterSequence.hh"
#include "fastjet/internal/MinHeap.hh"
#include "fastjet/internal/TilingExtent.hh"

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

using namespace std;


//----------------------------------------------------------------------
void ClusterSequence::_bj_remove_from_tiles(TiledJet * const jet) {
  Tile * tile = & _tiles[jet->tile_index];

  if (jet->previous == NULL) {
    // we are at head of the tile, so reset it.
    // If this was the only jet on the tile then tile->head will now be NULL
    tile->head = jet->next;
  } else {
    // adjust link from previous jet in this tile
    jet->previous->next = jet->next;
  }
  if (jet->next != NULL) {
    // adjust backwards-link from next jet in this tile
    jet->next->previous = jet->previous;
  }
}

//----------------------------------------------------------------------
/// Set up the tiles:
///  - decide the range in eta
///  - allocate the tiles
///  - set up the cross-referencing info between tiles
///
/// The neighbourhood of a tile is set up as follows
///
///           LRR
///           LXR
///           LLR
///
/// such that tiles is an array containing XLLLLRRRR with pointers
///                                         |   \ RH_tiles
///                                         \ surrounding_tiles
///
/// with appropriate precautions when close to the edge of the tiled
/// region.
///
void ClusterSequence::_initialise_tiles() {

  // first decide tile sizes (with a lower bound to avoid huge memory use with
  // very small R)
  double default_size = max(0.1,_Rparam);
  _tile_size_eta = default_size;
  // it makes no sense to go below 3 tiles in phi -- 3 tiles is
  // sufficient to make sure all pair-wise combinations up to pi in
  // phi are possible
  _n_tiles_phi   = max(3,int(floor(twopi/default_size)));
  _tile_size_phi = twopi / _n_tiles_phi; // >= _Rparam and fits in 2pi

  TilingExtent tiling_analysis(*this);
  _tiles_eta_min = tiling_analysis.minrap();
  _tiles_eta_max = tiling_analysis.maxrap();

  // // always include zero rapidity in the tiling region
  // _tiles_eta_min = 0.0;
  // _tiles_eta_max = 0.0;
  // // but go no further than following
  // const double maxrap = 7.0;
  // 
  // // and find out how much further one should go
  // for(unsigned int i = 0; i < _jets.size(); i++) {
  //   double eta = _jets[i].rap();
  //   // first check if eta is in range -- to avoid taking into account
  //   // very spurious rapidities due to particles with near-zero kt.
  //   if (abs(eta) < maxrap) {
  //     if (eta < _tiles_eta_min) {_tiles_eta_min = eta;}
  //     if (eta > _tiles_eta_max) {_tiles_eta_max = eta;}
  //   }
  // }

  // now adjust the values
  _tiles_ieta_min = int(floor(_tiles_eta_min/_tile_size_eta));
  _tiles_ieta_max = int(floor( _tiles_eta_max/_tile_size_eta));
  _tiles_eta_min = _tiles_ieta_min * _tile_size_eta;
  _tiles_eta_max = _tiles_ieta_max * _tile_size_eta;

  // allocate the tiles
  _tiles.resize((_tiles_ieta_max-_tiles_ieta_min+1)*_n_tiles_phi);

  // now set up the cross-referencing between tiles
  for (int ieta = _tiles_ieta_min; ieta <= _tiles_ieta_max; ieta++) {
    for (int iphi = 0; iphi < _n_tiles_phi; iphi++) {
      Tile * tile = & _tiles[_tile_index(ieta,iphi)];
      // no jets in this tile yet
      tile->head = NULL; // first element of tiles points to itself
      tile->begin_tiles[0] =  tile;
      Tile ** pptile = & (tile->begin_tiles[0]);
      pptile++;
      //
      // set up L's in column to the left of X
      tile->surrounding_tiles = pptile;
      if (ieta > _tiles_ieta_min) {
        // with the itile subroutine, we can safely run tiles from
        // idphi=-1 to idphi=+1, because it takes care of
        // negative and positive boundaries
        for (int idphi = -1; idphi <=+1; idphi++) {
          *pptile = & _tiles[_tile_index(ieta-1,iphi+idphi)];
          pptile++;
        }       
      }
      // now set up last L (below X)
      *pptile = & _tiles[_tile_index(ieta,iphi-1)];
      pptile++;
      // set up first R (above X)
      tile->RH_tiles = pptile;
      *pptile = & _tiles[_tile_index(ieta,iphi+1)];
      pptile++;
      // set up remaining R's, to the right of X
      if (ieta < _tiles_ieta_max) {
        for (int idphi = -1; idphi <= +1; idphi++) {
          *pptile = & _tiles[_tile_index(ieta+1,iphi+idphi)];
          pptile++;
        }       
      }
      // now put semaphore for end tile
      tile->end_tiles = pptile;
      // finally make sure tiles are untagged
      tile->tagged = false;
    }
  }

}


//----------------------------------------------------------------------
/// return the tile index corresponding to the given eta,phi point
int ClusterSequence::_tile_index(const double eta, const double phi) const {
  int ieta, iphi;
  if      (eta <= _tiles_eta_min) {ieta = 0;}
  else if (eta >= _tiles_eta_max) {ieta = _tiles_ieta_max-_tiles_ieta_min;}
  else {
    //ieta = int(floor((eta - _tiles_eta_min) / _tile_size_eta));
    ieta = int(((eta - _tiles_eta_min) / _tile_size_eta));
    // following needed in case of rare but nasty rounding errors
    if (ieta > _tiles_ieta_max-_tiles_ieta_min) {
      ieta = _tiles_ieta_max-_tiles_ieta_min;} 
  }
  // allow for some extent of being beyond range in calculation of phi
  // as well
  //iphi = (int(floor(phi/_tile_size_phi)) + _n_tiles_phi) % _n_tiles_phi;
  // with just int and no floor, things run faster but beware
  iphi = int((phi+twopi)/_tile_size_phi) % _n_tiles_phi;
  return (iphi + ieta * _n_tiles_phi);
}


//----------------------------------------------------------------------
// overloaded version which additionally sets up information regarding the
// tiling
inline void ClusterSequence::_tj_set_jetinfo( TiledJet * const jet,
                                              const int _jets_index) {
  // first call the generic setup
  _bj_set_jetinfo<>(jet, _jets_index);

  // Then do the setup specific to the tiled case.

  // Find out which tile it belonds to
  jet->tile_index = _tile_index(jet->eta, jet->phi);

  // Insert it into the tile's linked list of jets
  Tile * tile = &_tiles[jet->tile_index];
  jet->previous   = NULL;
  jet->next       = tile->head;
  if (jet->next != NULL) {jet->next->previous = jet;}
  tile->head      = jet;
}


//----------------------------------------------------------------------
/// output the contents of the tiles
void ClusterSequence::_print_tiles(TiledJet * briefjets ) const {
  for (vector<Tile>::const_iterator tile = _tiles.begin(); 
       tile < _tiles.end(); tile++) {
    cout << "Tile " << tile - _tiles.begin()<<" = ";
    vector<int> list;
    for (TiledJet * jetI = tile->head; jetI != NULL; jetI = jetI->next) {
      list.push_back(jetI-briefjets);
      //cout <<" "<<jetI-briefjets;
    }
    sort(list.begin(),list.end());
    for (unsigned int i = 0; i < list.size(); i++) {cout <<" "<<list[i];}
    cout <<"\n";
  }
}


//----------------------------------------------------------------------
/// Add to the vector tile_union the tiles that are in the neighbourhood
/// of the specified tile_index, including itself -- start adding
/// from position n_near_tiles-1, and increase n_near_tiles as
/// you go along (could have done it more C++ like with vector with reserved
/// space, but fear is that it would have been slower, e.g. checking
/// for end of vector at each stage to decide whether to resize it)
void ClusterSequence::_add_neighbours_to_tile_union(const int tile_index, 
               vector<int> & tile_union, int & n_near_tiles) const {
  for (Tile * const * near_tile = _tiles[tile_index].begin_tiles; 
       near_tile != _tiles[tile_index].end_tiles; near_tile++){
    // get the tile number
    tile_union[n_near_tiles] = *near_tile - & _tiles[0];
    n_near_tiles++;
  }
}


//----------------------------------------------------------------------
/// Like _add_neighbours_to_tile_union, but only adds neighbours if 
/// their "tagged" status is false; when a neighbour is added its
/// tagged status is set to true.
///
/// Note that with a high level of warnings (-pedantic -Wextra -ansi,
/// gcc complains about tile_index maybe being used uninitialised for
/// oldB in ClusterSequence::_minheap_faster_tiled_N2_cluster(). We
/// have explicitly checked that it was harmless so we could disable
/// the gcc warning by hand using the construct below
///
///  #pragma GCC diagnostic push
///  #pragma GCC diagnostic ignored "-Wpragmas"
///  #pragma GCC diagnostic ignored "-Wuninitialized"
///  #pragma GCC diagnostic ignored "-Wmaybe-uninitialized"
///    ...
///  #pragma GCC diagnostic pop
///
/// the @GCC diagnostic push/pop directive was only introduced in
/// gcc-4.6, so for broader usage, we'd need to insert #pragma GCC
/// diagnostic ignored "-Wpragmas" at the top of this file
inline void ClusterSequence::_add_untagged_neighbours_to_tile_union(
               const int tile_index, 
               vector<int> & tile_union, int & n_near_tiles)  {
  for (Tile ** near_tile = _tiles[tile_index].begin_tiles; 
       near_tile != _tiles[tile_index].end_tiles; near_tile++){
    if (! (*near_tile)->tagged) {
      (*near_tile)->tagged = true;
      // get the tile number
      tile_union[n_near_tiles] = *near_tile - & _tiles[0];
      n_near_tiles++;
    }
  }
}


//----------------------------------------------------------------------
/// run a tiled clustering
void ClusterSequence::_tiled_N2_cluster() {

  _initialise_tiles();

  int n = _jets.size();
  TiledJet * briefjets = new TiledJet[n];
  TiledJet * jetA = briefjets, * jetB;
  TiledJet oldB;
  oldB.tile_index=0; // prevents a gcc warning

  // will be used quite deep inside loops, but declare it here so that
  // memory (de)allocation gets done only once
  vector<int> tile_union(3*n_tile_neighbours);
  
  // initialise the basic jet info 
  for (int i = 0; i< n; i++) {
    _tj_set_jetinfo(jetA, i);
    //cout << i<<": "<<jetA->tile_index<<"\n";
    jetA++; // move on to next entry of briefjets
  }
  TiledJet * tail = jetA; // a semaphore for the end of briefjets
  TiledJet * head = briefjets; // a nicer way of naming start

  // set up the initial nearest neighbour information
  vector<Tile>::const_iterator tile;
  for (tile = _tiles.begin(); tile != _tiles.end(); tile++) {
    // first do it on this tile
    for (jetA = tile->head; jetA != NULL; jetA = jetA->next) {
      for (jetB = tile->head; jetB != jetA; jetB = jetB->next) {
        double dist = _bj_dist(jetA,jetB);
        if (dist < jetA->NN_dist) {jetA->NN_dist = dist; jetA->NN = jetB;}
        if (dist < jetB->NN_dist) {jetB->NN_dist = dist; jetB->NN = jetA;}
      }
    }
    // then do it for RH tiles
    for (Tile ** RTile = tile->RH_tiles; RTile != tile->end_tiles; RTile++) {
      for (jetA = tile->head; jetA != NULL; jetA = jetA->next) {
        for (jetB = (*RTile)->head; jetB != NULL; jetB = jetB->next) {
          double dist = _bj_dist(jetA,jetB);
          if (dist < jetA->NN_dist) {jetA->NN_dist = dist; jetA->NN = jetB;}
          if (dist < jetB->NN_dist) {jetB->NN_dist = dist; jetB->NN = jetA;}
        }
      }
    }
  }
  
  // now create the diJ (where J is i's NN) table -- remember that 
  // we differ from standard normalisation here by a factor of R2
  double * diJ = new double[n];
  jetA = head;
  for (int i = 0; i < n; i++) {
    diJ[i] = _bj_diJ(jetA);
    jetA++; // have jetA follow i
  }

  // now run the recombination loop
  int history_location = n-1;
  while (tail != head) {

    // find the minimum of the diJ on this round
    double diJ_min = diJ[0];
    int diJ_min_jet = 0;
    for (int i = 1; i < n; i++) {
      if (diJ[i] < diJ_min) {diJ_min_jet = i; diJ_min  = diJ[i];}
    }

    // do the recombination between A and B
    history_location++;
    jetA = & briefjets[diJ_min_jet];
    jetB = jetA->NN;
    // put the normalisation back in
    diJ_min *= _invR2; 

    //if (n == 19) {cout << "Hello "<<jetA-head<<" "<<jetB-head<<"\n";}

    //cout <<" WILL RECOMBINE "<< jetA-briefjets<<" "<<jetB-briefjets<<"\n";

    if (jetB != NULL) {
      // jet-jet recombination
      // If necessary relabel A & B to ensure jetB < jetA, that way if
      // the larger of them == newtail then that ends up being jetA and 
      // the new jet that is added as jetB is inserted in a position that
      // has a future!
      if (jetA < jetB) {std::swap(jetA,jetB);}

      int nn; // new jet index
      _do_ij_recombination_step(jetA->_jets_index, jetB->_jets_index, diJ_min, nn);

      //OBS// get the two history indices
      //OBSint hist_a = _jets[jetA->_jets_index].cluster_hist_index();
      //OBSint hist_b = _jets[jetB->_jets_index].cluster_hist_index();
      //OBS// create the recombined jet
      //OBS_jets.push_back(_jets[jetA->_jets_index] + _jets[jetB->_jets_index]);
      //OBSint nn = _jets.size() - 1;
      //OBS_jets[nn].set_cluster_hist_index(history_location);
      //OBS// update history
      //OBS//cout <<n-1<<" "<<jetA-head<<" "<<jetB-head<<"; ";
      //OBS_add_step_to_history(history_location, 
      //OBS                min(hist_a,hist_b),max(hist_a,hist_b),
      //OBS                        nn, diJ_min);

      // what was jetB will now become the new jet
      _bj_remove_from_tiles(jetA);
      oldB = * jetB;  // take a copy because we will need it...
      _bj_remove_from_tiles(jetB);
      _tj_set_jetinfo(jetB, nn); // also registers the jet in the tiling
    } else {
      // jet-beam recombination
      _do_iB_recombination_step(jetA->_jets_index, diJ_min);
            
      //OBS// get the hist_index
      //OBSint hist_a = _jets[jetA->_jets_index].cluster_hist_index();
      //OBS//cout <<n-1<<" "<<jetA-head<<" "<<-1<<"; ";
      //OBS_add_step_to_history(history_location,hist_a,BeamJet,Invalid,diJ_min); 
      _bj_remove_from_tiles(jetA);
    }

    // first establish the set of tiles over which we are going to 
    // have to run searches for updated and new nearest-neighbours
    int n_near_tiles = 0;
    _add_neighbours_to_tile_union(jetA->tile_index, tile_union, n_near_tiles);
    if (jetB != NULL) {
      bool sort_it = false;
      if (jetB->tile_index != jetA->tile_index) {
        sort_it = true;
        _add_neighbours_to_tile_union(jetB->tile_index,tile_union,n_near_tiles);
      }
      if (oldB.tile_index != jetA->tile_index && 
          oldB.tile_index != jetB->tile_index) {
        sort_it = true;
        _add_neighbours_to_tile_union(oldB.tile_index,tile_union,n_near_tiles);
      }

      if (sort_it) {
        // sort the tiles before then compressing the list
        sort(tile_union.begin(), tile_union.begin()+n_near_tiles);
        // and now condense the list
        int nnn = 1;
        for (int i = 1; i < n_near_tiles; i++) {
          if (tile_union[i] != tile_union[nnn-1]) {
            tile_union[nnn] = tile_union[i]; 
            nnn++;
          }
        }
        n_near_tiles = nnn;
      }
    }

    // now update our nearest neighbour info and diJ table
    // first reduce size of table
    tail--; n--;
    if (jetA == tail) {
      // there is nothing to be done
    } else {
      // Copy last jet contents and diJ info into position of jetA
      *jetA = *tail;
      diJ[jetA - head] = diJ[tail-head];
      // IN the tiling fix pointers to tail and turn them into
      // pointers to jetA (from predecessors, successors and the tile
      // head if need be)
      if (jetA->previous == NULL) {
        _tiles[jetA->tile_index].head = jetA;
      } else {
        jetA->previous->next = jetA;
      }
      if (jetA->next != NULL) {jetA->next->previous = jetA;}
    }

    // Initialise jetB's NN distance as well as updating it for 
    // other particles.
    for (int itile = 0; itile < n_near_tiles; itile++) {
      Tile * tile_ptr = &_tiles[tile_union[itile]];
      for (TiledJet * jetI = tile_ptr->head; jetI != NULL; jetI = jetI->next) {
        // see if jetI had jetA or jetB as a NN -- if so recalculate the NN
        if (jetI->NN == jetA || (jetI->NN == jetB && jetB != NULL)) {
          jetI->NN_dist = _R2;
          jetI->NN      = NULL;
          // now go over tiles that are neighbours of I (include own tile)
          for (Tile ** near_tile  = tile_ptr->begin_tiles; 
                       near_tile != tile_ptr->end_tiles; near_tile++) {
            // and then over the contents of that tile
            for (TiledJet * jetJ  = (*near_tile)->head; 
                            jetJ != NULL; jetJ = jetJ->next) {
              double dist = _bj_dist(jetI,jetJ);
              if (dist < jetI->NN_dist && jetJ != jetI) {
                jetI->NN_dist = dist; jetI->NN = jetJ;
              }
            }
          }
          diJ[jetI-head] = _bj_diJ(jetI); // update diJ 
        }
        // check whether new jetB is closer than jetI's current NN and
        // if need to update things
        if (jetB != NULL) {
          double dist = _bj_dist(jetI,jetB);
          if (dist < jetI->NN_dist) {
            if (jetI != jetB) {
              jetI->NN_dist = dist;
              jetI->NN = jetB;
              diJ[jetI-head] = _bj_diJ(jetI); // update diJ...
            }
          }
          if (dist < jetB->NN_dist) {
            if (jetI != jetB) {
              jetB->NN_dist = dist;
              jetB->NN      = jetI;}
          }
        }
      }
    }


    if (jetB != NULL) {diJ[jetB-head] = _bj_diJ(jetB);}
    //cout << n<<" "<<briefjets[95].NN-briefjets<<" "<<briefjets[95].NN_dist <<"\n";

    // remember to update pointers to tail
    for (Tile ** near_tile = _tiles[tail->tile_index].begin_tiles; 
                 near_tile!= _tiles[tail->tile_index].end_tiles; near_tile++){
      // and then the contents of that tile
      for (TiledJet * jetJ = (*near_tile)->head; 
                     jetJ != NULL; jetJ = jetJ->next) {
        if (jetJ->NN == tail) {jetJ->NN = jetA;}
      }
    }

    //for (int i = 0; i < n; i++) {
    //  if (briefjets[i].NN-briefjets >= n && briefjets[i].NN != NULL) {cout <<"YOU MUST BE CRAZY for n ="<<n<<", i = "<<i<<", NN = "<<briefjets[i].NN-briefjets<<"\n";}
    //}


    if (jetB != NULL) {diJ[jetB-head] = _bj_diJ(jetB);}
    //cout << briefjets[95].NN-briefjets<<" "<<briefjets[95].NN_dist <<"\n";

  }

  // final cleaning up;
  delete[] diJ;
  delete[] briefjets;
}


//----------------------------------------------------------------------
/// run a tiled clustering
void ClusterSequence::_faster_tiled_N2_cluster() {

  _initialise_tiles();

  int n = _jets.size();
  TiledJet * briefjets = new TiledJet[n];
  TiledJet * jetA = briefjets, * jetB;
  TiledJet oldB;
  oldB.tile_index=0; // prevents a gcc warning

  // will be used quite deep inside loops, but declare it here so that
  // memory (de)allocation gets done only once
  vector<int> tile_union(3*n_tile_neighbours);
  
  // initialise the basic jet info 
  for (int i = 0; i< n; i++) {
    _tj_set_jetinfo(jetA, i);
    //cout << i<<": "<<jetA->tile_index<<"\n";
    jetA++; // move on to next entry of briefjets
  }
  TiledJet * head = briefjets; // a nicer way of naming start

  // set up the initial nearest neighbour information
  vector<Tile>::const_iterator tile;
  for (tile = _tiles.begin(); tile != _tiles.end(); tile++) {
    // first do it on this tile
    for (jetA = tile->head; jetA != NULL; jetA = jetA->next) {
      for (jetB = tile->head; jetB != jetA; jetB = jetB->next) {
        double dist = _bj_dist(jetA,jetB);
        if (dist < jetA->NN_dist) {jetA->NN_dist = dist; jetA->NN = jetB;}
        if (dist < jetB->NN_dist) {jetB->NN_dist = dist; jetB->NN = jetA;}
      }
    }
    // then do it for RH tiles
    for (Tile ** RTile = tile->RH_tiles; RTile != tile->end_tiles; RTile++) {
      for (jetA = tile->head; jetA != NULL; jetA = jetA->next) {
        for (jetB = (*RTile)->head; jetB != NULL; jetB = jetB->next) {
          double dist = _bj_dist(jetA,jetB);
          if (dist < jetA->NN_dist) {jetA->NN_dist = dist; jetA->NN = jetB;}
          if (dist < jetB->NN_dist) {jetB->NN_dist = dist; jetB->NN = jetA;}
        }
      }
    }
    // no need to do it for LH tiles, since they are implicitly done
    // when we set NN for both jetA and jetB on the RH tiles.
  }
  
  // now create the diJ (where J is i's NN) table -- remember that 
  // we differ from standard normalisation here by a factor of R2
  // (corrected for at the end). 
  struct diJ_plus_link {
    double     diJ; // the distance
    TiledJet * jet; // the jet (i) for which we've found this distance
                    // (whose NN will the J).
  };
  diJ_plus_link * diJ = new diJ_plus_link[n];
  jetA = head;
  for (int i = 0; i < n; i++) {
    diJ[i].diJ = _bj_diJ(jetA); // kt distance * R^2
    diJ[i].jet = jetA;  // our compact diJ table will not be in      
    jetA->diJ_posn = i; // one-to-one corresp. with non-compact jets,
                        // so set up bi-directional correspondence here.
    jetA++; // have jetA follow i 
  }

  // now run the recombination loop
  int history_location = n-1;
  while (n > 0) {

    // find the minimum of the diJ on this round
    diJ_plus_link * best, *stop; // pointers a bit faster than indices
    // could use best to keep track of diJ min, but it turns out to be
    // marginally faster to have a separate variable (avoids n
    // dereferences at the expense of n/2 assignments).
    double diJ_min = diJ[0].diJ; // initialise the best one here.
    best = diJ;                  // and here
    stop = diJ+n;
    for (diJ_plus_link * here = diJ+1; here != stop; here++) {
      if (here->diJ < diJ_min) {best = here; diJ_min  = here->diJ;}
    }

    // do the recombination between A and B
    history_location++;
    jetA = best->jet;
    jetB = jetA->NN;
    // put the normalisation back in
    diJ_min *= _invR2; 

    if (jetB != NULL) {
      // jet-jet recombination
      // If necessary relabel A & B to ensure jetB < jetA, that way if
      // the larger of them == newtail then that ends up being jetA and 
      // the new jet that is added as jetB is inserted in a position that
      // has a future!
      if (jetA < jetB) {std::swap(jetA,jetB);}

      int nn; // new jet index
      _do_ij_recombination_step(jetA->_jets_index, jetB->_jets_index, diJ_min, nn);
      
      //OBS// get the two history indices
      //OBSint ihstry_a = _jets[jetA->_jets_index].cluster_hist_index();
      //OBSint ihstry_b = _jets[jetB->_jets_index].cluster_hist_index();
      //OBS// create the recombined jet
      //OBS_jets.push_back(_jets[jetA->_jets_index] + _jets[jetB->_jets_index]);
      //OBSint nn = _jets.size() - 1;
      //OBS_jets[nn].set_cluster_hist_index(history_location);
      //OBS// update history
      //OBS//cout <<n-1<<" "<<jetA-head<<" "<<jetB-head<<"; ";
      //OBS_add_step_to_history(history_location, 
      //OBS                min(ihstry_a,ihstry_b),max(ihstry_a,ihstry_b),
      //OBS                        nn, diJ_min);
      // what was jetB will now become the new jet
      _bj_remove_from_tiles(jetA);
      oldB = * jetB;  // take a copy because we will need it...
      _bj_remove_from_tiles(jetB);
      _tj_set_jetinfo(jetB, nn); // cause jetB to become _jets[nn]
                                 // (also registers the jet in the tiling)
    } else {
      // jet-beam recombination
      // get the hist_index
      _do_iB_recombination_step(jetA->_jets_index, diJ_min);
      //OBSint ihstry_a = _jets[jetA->_jets_index].cluster_hist_index();
      //OBS//cout <<n-1<<" "<<jetA-head<<" "<<-1<<"; ";
      //OBS_add_step_to_history(history_location,ihstry_a,BeamJet,Invalid,diJ_min); 
      _bj_remove_from_tiles(jetA);
    }

    // first establish the set of tiles over which we are going to
    // have to run searches for updated and new nearest-neighbours --
    // basically a combination of vicinity of the tiles of the two old
    // and one new jet.
    int n_near_tiles = 0;
    _add_untagged_neighbours_to_tile_union(jetA->tile_index, 
                                           tile_union, n_near_tiles);
    if (jetB != NULL) {
      if (jetB->tile_index != jetA->tile_index) {
        _add_untagged_neighbours_to_tile_union(jetB->tile_index,
                                               tile_union,n_near_tiles);
      }
      if (oldB.tile_index != jetA->tile_index && 
          oldB.tile_index != jetB->tile_index) {
        _add_untagged_neighbours_to_tile_union(oldB.tile_index,
                                               tile_union,n_near_tiles);
      }
    }

    // now update our nearest neighbour info and diJ table
    // first reduce size of table
    n--;
    // then compactify the diJ by taking the last of the diJ and copying
    // it to the position occupied by the diJ for jetA
    diJ[n].jet->diJ_posn = jetA->diJ_posn;
    diJ[jetA->diJ_posn] = diJ[n];

    // Initialise jetB's NN distance as well as updating it for 
    // other particles.
    // Run over all tiles in our union 
    for (int itile = 0; itile < n_near_tiles; itile++) {
      Tile * tile_ptr = &_tiles[tile_union[itile]];
      tile_ptr->tagged = false; // reset tag, since we're done with unions
      // run over all jets in the current tile
      for (TiledJet * jetI = tile_ptr->head; jetI != NULL; jetI = jetI->next) {
        // see if jetI had jetA or jetB as a NN -- if so recalculate the NN
        if (jetI->NN == jetA || (jetI->NN == jetB && jetB != NULL)) {
          jetI->NN_dist = _R2;
          jetI->NN      = NULL;
          // now go over tiles that are neighbours of I (include own tile)
          for (Tile ** near_tile  = tile_ptr->begin_tiles; 
                       near_tile != tile_ptr->end_tiles; near_tile++) {
            // and then over the contents of that tile
            for (TiledJet * jetJ  = (*near_tile)->head; 
                            jetJ != NULL; jetJ = jetJ->next) {
              double dist = _bj_dist(jetI,jetJ);
              if (dist < jetI->NN_dist && jetJ != jetI) {
                jetI->NN_dist = dist; jetI->NN = jetJ;
              }
            }
          }
          diJ[jetI->diJ_posn].diJ = _bj_diJ(jetI); // update diJ kt-dist
        }
        // check whether new jetB is closer than jetI's current NN and
        // if jetI is closer than jetB's current (evolving) nearest
        // neighbour. Where relevant update things
        if (jetB != NULL) {
          double dist = _bj_dist(jetI,jetB);
          if (dist < jetI->NN_dist) {
            if (jetI != jetB) {
              jetI->NN_dist = dist;
              jetI->NN = jetB;
              diJ[jetI->diJ_posn].diJ = _bj_diJ(jetI); // update diJ...
            }
          }
          if (dist < jetB->NN_dist) {
            if (jetI != jetB) {
              jetB->NN_dist = dist;
              jetB->NN      = jetI;}
          }
        }
      }
    }

    // finally, register the updated kt distance for B
    if (jetB != NULL) {diJ[jetB->diJ_posn].diJ = _bj_diJ(jetB);}

  }

  // final cleaning up;
  delete[] diJ;
  delete[] briefjets;
}

//----------------------------------------------------------------------
/// run a tiled clustering, with our minheap for keeping track of the
/// smallest dij
void ClusterSequence::_minheap_faster_tiled_N2_cluster() {

  _initialise_tiles();

  int n = _jets.size();
  TiledJet * briefjets = new TiledJet[n];
  TiledJet * jetA = briefjets, * jetB;
  TiledJet oldB;
  oldB.tile_index=0; // prevents a gcc warning

  // will be used quite deep inside loops, but declare it here so that
  // memory (de)allocation gets done only once
  vector<int> tile_union(3*n_tile_neighbours);
  
  // initialise the basic jet info 
  for (int i = 0; i< n; i++) {
    _tj_set_jetinfo(jetA, i);
    //cout << i<<": "<<jetA->tile_index<<"\n";
    jetA++; // move on to next entry of briefjets
  }
  TiledJet * head = briefjets; // a nicer way of naming start

  // set up the initial nearest neighbour information
  vector<Tile>::const_iterator tile;
  for (tile = _tiles.begin(); tile != _tiles.end(); tile++) {
    // first do it on this tile
    for (jetA = tile->head; jetA != NULL; jetA = jetA->next) {
      for (jetB = tile->head; jetB != jetA; jetB = jetB->next) {
        double dist = _bj_dist(jetA,jetB);
        if (dist < jetA->NN_dist) {jetA->NN_dist = dist; jetA->NN = jetB;}
        if (dist < jetB->NN_dist) {jetB->NN_dist = dist; jetB->NN = jetA;}
      }
    }
    // then do it for RH tiles
    for (Tile ** RTile = tile->RH_tiles; RTile != tile->end_tiles; RTile++) {
      for (jetA = tile->head; jetA != NULL; jetA = jetA->next) {
        for (jetB = (*RTile)->head; jetB != NULL; jetB = jetB->next) {
          double dist = _bj_dist(jetA,jetB);
          if (dist < jetA->NN_dist) {jetA->NN_dist = dist; jetA->NN = jetB;}
          if (dist < jetB->NN_dist) {jetB->NN_dist = dist; jetB->NN = jetA;}
        }
      }
    }
    // no need to do it for LH tiles, since they are implicitly done
    // when we set NN for both jetA and jetB on the RH tiles.
  }

  
  //// now create the diJ (where J is i's NN) table -- remember that 
  //// we differ from standard normalisation here by a factor of R2
  //// (corrected for at the end). 
  //struct diJ_plus_link {
  //  double     diJ; // the distance
  //  TiledJet * jet; // the jet (i) for which we've found this distance
  //                  // (whose NN will the J).
  //};
  //diJ_plus_link * diJ = new diJ_plus_link[n];
  //jetA = head;
  //for (int i = 0; i < n; i++) {
  //  diJ[i].diJ = _bj_diJ(jetA); // kt distance * R^2
  //  diJ[i].jet = jetA;  // our compact diJ table will not be in            
  //  jetA->diJ_posn = i; // one-to-one corresp. with non-compact jets,
  //                      // so set up bi-directional correspondence here.
  //  jetA++; // have jetA follow i 
  //}

  vector<double> diJs(n);
  for (int i = 0; i < n; i++) {
    diJs[i] = _bj_diJ(&briefjets[i]);
    briefjets[i].label_minheap_update_done();
  }
  MinHeap minheap(diJs);
  // have a stack telling us which jets we'll have to update on the heap
  vector<TiledJet *> jets_for_minheap;
  jets_for_minheap.reserve(n); 

  // now run the recombination loop
  int history_location = n-1;
  while (n > 0) {

    double diJ_min = minheap.minval() *_invR2;
    jetA = head + minheap.minloc();

    // do the recombination between A and B
    history_location++;
    jetB = jetA->NN;

    if (jetB != NULL) {
      // jet-jet recombination
      // If necessary relabel A & B to ensure jetB < jetA, that way if
      // the larger of them == newtail then that ends up being jetA and 
      // the new jet that is added as jetB is inserted in a position that
      // has a future!
      if (jetA < jetB) {std::swap(jetA,jetB);}

      int nn; // new jet index
      _do_ij_recombination_step(jetA->_jets_index, jetB->_jets_index, diJ_min, nn);
      
      // what was jetB will now become the new jet
      _bj_remove_from_tiles(jetA);
      oldB = * jetB;  // take a copy because we will need it...
      _bj_remove_from_tiles(jetB);
      _tj_set_jetinfo(jetB, nn); // cause jetB to become _jets[nn]
                                 // (also registers the jet in the tiling)
    } else {
      // jet-beam recombination
      // get the hist_index
      _do_iB_recombination_step(jetA->_jets_index, diJ_min);
      _bj_remove_from_tiles(jetA);
    }

    // remove the minheap entry for jetA
    minheap.remove(jetA-head);

    // first establish the set of tiles over which we are going to
    // have to run searches for updated and new nearest-neighbours --
    // basically a combination of vicinity of the tiles of the two old
    // and one new jet.
    int n_near_tiles = 0;
    _add_untagged_neighbours_to_tile_union(jetA->tile_index, 
                                           tile_union, n_near_tiles);
    if (jetB != NULL) {
      if (jetB->tile_index != jetA->tile_index) {
        _add_untagged_neighbours_to_tile_union(jetB->tile_index,
                                               tile_union,n_near_tiles);
      }
      if (oldB.tile_index != jetA->tile_index && 
          oldB.tile_index != jetB->tile_index) {
        // GS: the line below generates a warning that oldB.tile_index
        // may be used uninitialised. However, to reach this point, we
        // ned jetB != NULL (see test a few lines above) and is jetB
        // !=NULL, one would have gone through "oldB = *jetB before
        // (see piece of code ~20 line above), so the index is
        // initialised. We do not do anything to avoid the warning to
        // avoid any potential speed impact.
        _add_untagged_neighbours_to_tile_union(oldB.tile_index,
                                               tile_union,n_near_tiles);
      }
      // indicate that we'll have to update jetB in the minheap
      jetB->label_minheap_update_needed();
      jets_for_minheap.push_back(jetB);
    }


    // Initialise jetB's NN distance as well as updating it for 
    // other particles.
    // Run over all tiles in our union 
    for (int itile = 0; itile < n_near_tiles; itile++) {
      Tile * tile_ptr = &_tiles[tile_union[itile]];
      tile_ptr->tagged = false; // reset tag, since we're done with unions
      // run over all jets in the current tile
      for (TiledJet * jetI = tile_ptr->head; jetI != NULL; jetI = jetI->next) {
        // see if jetI had jetA or jetB as a NN -- if so recalculate the NN
        if (jetI->NN == jetA || (jetI->NN == jetB && jetB != NULL)) {
          jetI->NN_dist = _R2;
          jetI->NN      = NULL;
          // label jetI as needing heap action...
          if (!jetI->minheap_update_needed()) {
            jetI->label_minheap_update_needed();
            jets_for_minheap.push_back(jetI);}
          // now go over tiles that are neighbours of I (include own tile)
          for (Tile ** near_tile  = tile_ptr->begin_tiles; 
                       near_tile != tile_ptr->end_tiles; near_tile++) {
            // and then over the contents of that tile
            for (TiledJet * jetJ  = (*near_tile)->head; 
                            jetJ != NULL; jetJ = jetJ->next) {
              double dist = _bj_dist(jetI,jetJ);
              if (dist < jetI->NN_dist && jetJ != jetI) {
                jetI->NN_dist = dist; jetI->NN = jetJ;
              }
            }
          }
        }
        // check whether new jetB is closer than jetI's current NN and
        // if jetI is closer than jetB's current (evolving) nearest
        // neighbour. Where relevant update things
        if (jetB != NULL) {
          double dist = _bj_dist(jetI,jetB);
          if (dist < jetI->NN_dist) {
            if (jetI != jetB) {
              jetI->NN_dist = dist;
              jetI->NN = jetB;
              // label jetI as needing heap action...
              if (!jetI->minheap_update_needed()) {
                jetI->label_minheap_update_needed();
                jets_for_minheap.push_back(jetI);}
            }
          }
          if (dist < jetB->NN_dist) {
            if (jetI != jetB) {
              jetB->NN_dist = dist;
              jetB->NN      = jetI;}
          }
        }
      }
    }

    // deal with jets whose minheap entry needs updating
    while (jets_for_minheap.size() > 0) {
      TiledJet * jetI = jets_for_minheap.back(); 
      jets_for_minheap.pop_back();
      minheap.update(jetI-head, _bj_diJ(jetI));
      jetI->label_minheap_update_done();
    }
    n--;
  }

  // final cleaning up;
  delete[] briefjets;
}


FASTJET_END_NAMESPACE