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