Loading...
Searching...
No Matches
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-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
//----------------------------------------------------------------------
/// 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 // __FASTJET_CLUSTERQUENCE_N2_ICC__