ClusterSequence_N2.icc

// -*- C++ -*-
#ifndef __FASTJET_CLUSTERQUENCE_N2_ICC__
#define __FASTJET_CLUSTERQUENCE_N2_ICC__
#include "fastjet/ClusterSequence.hh"
 
//FJSTARTHEADER
// $Id: ClusterSequence_N2.cc 1351 2009-01-09 18:03:03Z salam $
//
// Copyright (c) 2005-2025, 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
 
//----------------------------------------------------------------------
/// 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
  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
    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.
 
    // There are several potential variants of the algorithm below for
    // updating the NN info and diJ table. Differences affect cases
    // where jetB is null. Numbers below are for antikt, R=0.4,
    // Nparticles=50, cf. top-level/development-notes/2025-03-speed and
    // the percent number indicates the change in runtime relative to
    // "original" (more negative is better)
    //
    // 0: original, i.e. up to FJ version 3.4.3, good compromise in normal conditions, pathological for very small R (cf. 2503.08146)
    // 1: add extra non-null jetB test in loop             intel/linux:+1%, M2Pro/mac:+2%
    // 2: two independent loops                            intel/linux:-3%, M2Pro/mac:+2%
    // 3: proxy jetB_or_A to avoid null test               intel/linux:-1%, M2Pro/mac:-1%
    //
    // Note that on intel/linux/gcc(14.2) -O2 v. -O3 tends to be a small difference
    // while on M2Pro/mac/clang(17) -O3 is about 10% faster than -O2, while gcc is 7% faster and less sensitive to -O2 v. -O3
    // 
    // We settle on strategy 2 for the 3.5.X series, starting 2025-05, on the
    // assumption that intel/linux is the dominant platform for most
    // very long runs with FJ.
#define FJ_N2_UPDATE_STRATEGY 2
 
#if   FJ_N2_UPDATE_STRATEGY==0
    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);}
#elif FJ_N2_UPDATE_STRATEGY==1
    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 || (jetB && 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);}
#elif FJ_N2_UPDATE_STRATEGY==2
    if (jetB != NULL){
      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
        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...
          }
        }
        // and check whether jetI is potentially jetB's NN
        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;}
      }
      diJ[jetB-head] = _bj_diJ(jetB);
    } else {
      // case where jetB == NULL
      for (BJ * jetI = head; jetI != tail; jetI++) {
        // see if jetI had jetA as a NN -- if so recalculate the NN
        if (jetI->NN == jetA) {
          _bj_set_NN_nocross(jetI, head, tail);
          diJ[jetI-head] = _bj_diJ(jetI); // update diJ 
        } 
        // if jetI's NN is the new tail then relabel it so that it becomes jetA
        if (jetI->NN == tail) {jetI->NN = jetA;}
      }
    }
#elif FJ_N2_UPDATE_STRATEGY==3
    auto jetB_or_A = jetB ? jetB : jetA;
    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_or_A) {
        _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);}
#else
#error "unknown N2 update strategy"
#endif
#undef FJ_N2_UPDATE_STRATEGY
  }
  
 
  // 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 // __FASTJET_CLUSTERQUENCE_N2_ICC__