FastJet 3.0.4
ClusterSequence_N2.icc
00001 // -*- C++ -*-
00002 #ifndef __CLUSTERQUENCE_N2_ICC__
00003 #define __CLUSTERQUENCE_N2_ICC__
00004 #include "fastjet/ClusterSequence.hh"
00005 
00006 //STARTHEADER
00007 // $Id: ClusterSequence_N2.cc 1351 2009-01-09 18:03:03Z salam $
00008 //
00009 // Copyright (c) 2005-2009, Matteo Cacciari, Gavin Salam and Gregory Soyez
00010 //
00011 //----------------------------------------------------------------------
00012 // This file is part of FastJet.
00013 //
00014 //  FastJet is free software; you can redistribute it and/or modify
00015 //  it under the terms of the GNU General Public License as published by
00016 //  the Free Software Foundation; either version 2 of the License, or
00017 //  (at your option) any later version.
00018 //
00019 //  The algorithms that underlie FastJet have required considerable
00020 //  development and are described in hep-ph/0512210. If you use
00021 //  FastJet as part of work towards a scientific publication, please
00022 //  include a citation to the FastJet paper.
00023 //
00024 //  FastJet is distributed in the hope that it will be useful,
00025 //  but WITHOUT ANY WARRANTY; without even the implied warranty of
00026 //  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
00027 //  GNU General Public License for more details.
00028 //
00029 //  You should have received a copy of the GNU General Public License
00030 //  along with FastJet; if not, write to the Free Software
00031 //  Foundation, Inc.:
00032 //      59 Temple Place, Suite 330, Boston, MA  02111-1307  USA
00033 //----------------------------------------------------------------------
00034 //ENDHEADER
00035 
00036 //----------------------------------------------------------------------
00037 /// Order(N^2) clustering 
00038 ///
00039 /// Works for any class BJ that satisfies certain minimal 
00040 /// requirements (which are ...?)
00041 ///
00042 /// - need to have _bj_set_jetinfo
00043 /// - need to have _bj_dist
00044 /// - should contain members kt2 (=energy^2), NN, NN_dist, _jets_index
00045 
00046 FASTJET_BEGIN_NAMESPACE      // defined in fastjet/internal/base.hh
00047 
00048 // this should not normally appear in an include file, but we'll make an 
00049 // exception seeing as this is
00050 //using namespace std;
00051 
00052 
00053 template<class BJ> void ClusterSequence::_simple_N2_cluster() {
00054   int n = _jets.size();
00055   BJ * briefjets = new BJ[n];
00056   BJ * jetA = briefjets, * jetB;
00057   
00058   // initialise the basic jet info 
00059   for (int i = 0; i< n; i++) {
00060     _bj_set_jetinfo(jetA, i);
00061     jetA++; // move on to next entry of briefjets
00062   }
00063   BJ * tail = jetA; // a semaphore for the end of briefjets
00064   BJ * head = briefjets; // a nicer way of naming start
00065 
00066   // now initialise the NN distances: jetA will run from 1..n-1; and
00067   // jetB from 0..jetA-1
00068   for (jetA = head + 1; jetA != tail; jetA++) {
00069     // set NN info for jetA based on jets running from head..jetA-1,
00070     // checking in the process whether jetA itself is an undiscovered
00071     // NN of one of those jets.
00072     _bj_set_NN_crosscheck(jetA, head, jetA);
00073   }
00074 
00075 
00076   // now create the diJ (where J is i's NN) table -- remember that 
00077   // we differ from standard normalisation here by a factor of R2
00078   double * diJ = new double[n];
00079   jetA = head;
00080   for (int i = 0; i < n; i++) {
00081     diJ[i] = _bj_diJ(jetA);
00082     jetA++; // have jetA follow i
00083   }
00084 
00085   // now run the recombination loop
00086   int history_location = n-1;
00087   while (tail != head) {
00088 
00089     // find the minimum of the diJ on this round
00090     double diJ_min = diJ[0];
00091     int diJ_min_jet = 0;
00092     for (int i = 1; i < n; i++) {
00093       if (diJ[i] < diJ_min) {diJ_min_jet = i; diJ_min  = diJ[i];}
00094     }
00095     
00096     // do the recombination between A and B
00097     history_location++;
00098     jetA = & briefjets[diJ_min_jet];
00099     // GPS mod 2009-02-11
00100     //jetB = jetA->NN;
00101     jetB = static_cast<BJ *>(jetA->NN);
00102     // put the normalisation back in
00103     diJ_min *= _invR2; 
00104     if (jetB != NULL) {
00105       // jet-jet recombination
00106       // If necessary relabel A & B to ensure jetB < jetA, that way if
00107       // the larger of them == newtail then that ends up being jetA and 
00108       // the new jet that is added as jetB is inserted in a position that
00109       // has a future!
00110       if (jetA < jetB) {std::swap(jetA,jetB);}
00111 
00112       int nn; // new jet index
00113       _do_ij_recombination_step(jetA->_jets_index, jetB->_jets_index, diJ_min, nn);
00114 
00115       // what was jetB will now become the new jet
00116       _bj_set_jetinfo(jetB, nn);
00117 
00118     } else {
00119       // jet-beam recombination
00120       _do_iB_recombination_step(jetA->_jets_index, diJ_min);
00121     }
00122 
00123     // now update our nearest neighbour info and diJ table
00124     // first reduce size of table
00125     tail--; n--;
00126     // Copy last jet contents and diJ info into position of jetA
00127     *jetA = *tail;
00128     diJ[jetA - head] = diJ[tail-head];
00129 
00130     // Initialise jetB's NN distance as well as updating it for 
00131     // other particles.
00132     // NB: by having different loops for jetB == or != NULL we could
00133     //     perhaps save a few percent (usually avoid one if inside loop),
00134     //     but will not do it for now because on laptop fluctuations are
00135     //     too large to reliably measure a few percent difference...
00136     for (BJ * jetI = head; jetI != tail; jetI++) {
00137       // see if jetI had jetA or jetB as a NN -- if so recalculate the NN
00138       if (jetI->NN == jetA || jetI->NN == jetB) {
00139         _bj_set_NN_nocross(jetI, head, tail);
00140         diJ[jetI-head] = _bj_diJ(jetI); // update diJ 
00141       } 
00142       // check whether new jetB is closer than jetI's current NN and
00143       // if need be update things
00144       if (jetB != NULL) {
00145         double dist = _bj_dist(jetI,jetB);
00146         if (dist < jetI->NN_dist) {
00147           if (jetI != jetB) {
00148             jetI->NN_dist = dist;
00149             jetI->NN = jetB;
00150             diJ[jetI-head] = _bj_diJ(jetI); // update diJ...
00151           }
00152         }
00153         if (dist < jetB->NN_dist) {
00154           if (jetI != jetB) {
00155             jetB->NN_dist = dist;
00156             jetB->NN      = jetI;}
00157         }
00158       }
00159       // if jetI's NN is the new tail then relabel it so that it becomes jetA
00160       if (jetI->NN == tail) {jetI->NN = jetA;}
00161     }
00162 
00163     
00164     if (jetB != NULL) {diJ[jetB-head] = _bj_diJ(jetB);}
00165 
00166 
00167   }
00168   
00169 
00170   // final cleaning up;
00171   delete[] diJ;
00172   delete[] briefjets;
00173 }
00174 
00175 
00176 
00177 // //----------------------------------------------------------------------
00178 // // initialises a GenBriefJet
00179 // template<> inline void ClusterSequence::_bj_set_jetinfo(
00180 //                            GenBriefJet * const jetA, const int _jets_index) const {
00181 // 
00182 //   jetA->init(_jets[_jets_index]);
00183 //   jetA->_jets_index = _jets_index;
00184 // 
00185 // }
00186 // 
00187 // 
00188 // //----------------------------------------------------------------------
00189 // // returns the distance between two GenBriefJets
00190 // template<> double ClusterSequence::_bj_dist(
00191 //                 const GenBriefJet * const jeta, 
00192 //                 const GenBriefJet * const jetb) const {
00193 //   return jeta->geom_ij(jetb);
00194 // }
00195 
00196 FASTJET_END_NAMESPACE
00197 
00198 #endif // __CLUSTERQUENCE_N2_ICC__
00199 
00200 
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends