ClusterSequence_N2.icc

// -*- C++ -*-
#ifndef __CLUSTERQUENCE_N2_ICC__
#define __CLUSTERQUENCE_N2_ICC__
#include "fastjet/ClusterSequence.hh"

//STARTHEADER
// $Id: ClusterSequence_N2.cc 1351 2009-01-09 18:03:03Z salam $
//
// Copyright (c) 2005-2009, Matteo Cacciari, Gavin 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 and are described in hep-ph/0512210. If you use
//  FastJet as part of work towards a scientific publication, please
//  include a citation to the FastJet paper.
//
//  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, write to the Free Software
//  Foundation, Inc.:
//      59 Temple Place, Suite 330, Boston, MA  02111-1307  USA
//----------------------------------------------------------------------
//ENDHEADER

//----------------------------------------------------------------------
/// Order(N^2) clustering 
///
/// Works for any class BJ that satisfies certain minimal 
/// requirements (which are ...?)
///
/// - need to have _bj_set_jetinfo
/// - need to have _bj_dist
/// - should contain members kt2 (=energy^2), NN, NN_dist, _jets_index

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

// this should not normally appear in an include file, but we'll make an 
// exception seeing as this is
//using namespace std;


template<class BJ> void ClusterSequence::_simple_N2_cluster() {
  int n = _jets.size();
  BJ * briefjets = new BJ[n];
  BJ * jetA = briefjets, * jetB;
  
  // initialise the basic jet info 
  for (int i = 0; i< n; i++) {
    _bj_set_jetinfo(jetA, i);
    jetA++; // move on to next entry of briefjets
  }
  BJ * tail = jetA; // a semaphore for the end of briefjets
  BJ * head = briefjets; // a nicer way of naming start

  // now initialise the NN distances: jetA will run from 1..n-1; and
  // jetB from 0..jetA-1
  for (jetA = head + 1; jetA != tail; jetA++) {
    // set NN info for jetA based on jets running from head..jetA-1,
    // checking in the process whether jetA itself is an undiscovered
    // NN of one of those jets.
    _bj_set_NN_crosscheck(jetA, head, 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];
    // GPS mod 2009-02-11
    //jetB = jetA->NN;
    jetB = static_cast<BJ *>(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);

      // what was jetB will now become the new jet
      _bj_set_jetinfo(jetB, nn);

    } else {
      // jet-beam recombination
      _do_iB_recombination_step(jetA->_jets_index, diJ_min);
    }

    // now update our nearest neighbour info and diJ table
    // first reduce size of table
    tail--; n--;
    // Copy last jet contents and diJ info into position of jetA
    *jetA = *tail;
    diJ[jetA - head] = diJ[tail-head];

    // Initialise jetB's NN distance as well as updating it for 
    // other particles.
    // NB: by having different loops for jetB == or != NULL we could
    //     perhaps save a few percent (usually avoid one if inside loop),
    //     but will not do it for now because on laptop fluctuations are
    //     too large to reliably measure a few percent difference...
    for (BJ * jetI = head; jetI != tail; jetI++) {
      // see if jetI had jetA or jetB as a NN -- if so recalculate the NN
      if (jetI->NN == jetA || jetI->NN == jetB) {
        _bj_set_NN_nocross(jetI, head, tail);
        diJ[jetI-head] = _bj_diJ(jetI); // update diJ 
      } 
      // check whether new jetB is closer than jetI's current NN and
      // if need be 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 jetI's NN is the new tail then relabel it so that it becomes jetA
      if (jetI->NN == tail) {jetI->NN = jetA;}
    }

    
    if (jetB != NULL) {diJ[jetB-head] = _bj_diJ(jetB);}


  }
  

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



// //----------------------------------------------------------------------
// // initialises a GenBriefJet
// template<> inline void ClusterSequence::_bj_set_jetinfo(
//                            GenBriefJet * const jetA, const int _jets_index) const {
// 
//   jetA->init(_jets[_jets_index]);
//   jetA->_jets_index = _jets_index;
// 
// }
// 
// 
// //----------------------------------------------------------------------
// // returns the distance between two GenBriefJets
// template<> double ClusterSequence::_bj_dist(
//                 const GenBriefJet * const jeta, 
//                 const GenBriefJet * const jetb) const {
//   return jeta->geom_ij(jetb);
// }

FASTJET_END_NAMESPACE

#endif // __CLUSTERQUENCE_N2_ICC__