FastJet  3.1.3
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Modules Pages
ClusterSequence_N2.icc
1 // -*- C++ -*-
2 #ifndef __FASTJET_CLUSTERQUENCE_N2_ICC__
3 #define __FASTJET_CLUSTERQUENCE_N2_ICC__
4 #include "fastjet/ClusterSequence.hh"
5 
6 //FJSTARTHEADER
7 // $Id: ClusterSequence_N2.cc 1351 2009-01-09 18:03:03Z salam $
8 //
9 // Copyright (c) 2005-2014, Matteo Cacciari, Gavin P. Salam and Gregory Soyez
10 //
11 //----------------------------------------------------------------------
12 // This file is part of FastJet.
13 //
14 // FastJet is free software; you can redistribute it and/or modify
15 // it under the terms of the GNU General Public License as published by
16 // the Free Software Foundation; either version 2 of the License, or
17 // (at your option) any later version.
18 //
19 // The algorithms that underlie FastJet have required considerable
20 // development. They are described in the original FastJet paper,
21 // hep-ph/0512210 and in the manual, arXiv:1111.6097. If you use
22 // FastJet as part of work towards a scientific publication, please
23 // quote the version you use and include a citation to the manual and
24 // optionally also to hep-ph/0512210.
25 //
26 // FastJet is distributed in the hope that it will be useful,
27 // but WITHOUT ANY WARRANTY; without even the implied warranty of
28 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
29 // GNU General Public License for more details.
30 //
31 // You should have received a copy of the GNU General Public License
32 // along with FastJet. If not, see <http://www.gnu.org/licenses/>.
33 //----------------------------------------------------------------------
34 //FJENDHEADER
35 
36 //----------------------------------------------------------------------
37 /// Order(N^2) clustering
38 ///
39 /// Works for any class BJ that satisfies certain minimal
40 /// requirements (which are ...?)
41 ///
42 /// - need to have _bj_set_jetinfo
43 /// - need to have _bj_dist
44 /// - should contain members kt2 (=energy^2), NN, NN_dist, _jets_index
45 
46 FASTJET_BEGIN_NAMESPACE // defined in fastjet/internal/base.hh
47 
48 // this should not normally appear in an include file, but we'll make an
49 // exception seeing as this is
50 //using namespace std;
51 
52 
53 template<class BJ> void ClusterSequence::_simple_N2_cluster() {
54  int n = _jets.size();
55  BJ * briefjets = new BJ[n];
56  BJ * jetA = briefjets, * jetB;
57 
58  // initialise the basic jet info
59  for (int i = 0; i< n; i++) {
60  _bj_set_jetinfo(jetA, i);
61  jetA++; // move on to next entry of briefjets
62  }
63  BJ * tail = jetA; // a semaphore for the end of briefjets
64  BJ * head = briefjets; // a nicer way of naming start
65 
66  // now initialise the NN distances: jetA will run from 1..n-1; and
67  // jetB from 0..jetA-1
68  for (jetA = head + 1; jetA != tail; jetA++) {
69  // set NN info for jetA based on jets running from head..jetA-1,
70  // checking in the process whether jetA itself is an undiscovered
71  // NN of one of those jets.
72  _bj_set_NN_crosscheck(jetA, head, jetA);
73  }
74 
75 
76  // now create the diJ (where J is i's NN) table -- remember that
77  // we differ from standard normalisation here by a factor of R2
78  double * diJ = new double[n];
79  jetA = head;
80  for (int i = 0; i < n; i++) {
81  diJ[i] = _bj_diJ(jetA);
82  jetA++; // have jetA follow i
83  }
84 
85  // now run the recombination loop
86  int history_location = n-1;
87  while (tail != head) {
88 
89  // find the minimum of the diJ on this round
90  double diJ_min = diJ[0];
91  int diJ_min_jet = 0;
92  for (int i = 1; i < n; i++) {
93  if (diJ[i] < diJ_min) {diJ_min_jet = i; diJ_min = diJ[i];}
94  }
95 
96  // do the recombination between A and B
97  history_location++;
98  jetA = & briefjets[diJ_min_jet];
99  // GPS mod 2009-02-11
100  //jetB = jetA->NN;
101  jetB = static_cast<BJ *>(jetA->NN);
102  // put the normalisation back in
103  diJ_min *= _invR2;
104  if (jetB != NULL) {
105  // jet-jet recombination
106  // If necessary relabel A & B to ensure jetB < jetA, that way if
107  // the larger of them == newtail then that ends up being jetA and
108  // the new jet that is added as jetB is inserted in a position that
109  // has a future!
110  if (jetA < jetB) {std::swap(jetA,jetB);}
111 
112  int nn; // new jet index
113  _do_ij_recombination_step(jetA->_jets_index, jetB->_jets_index, diJ_min, nn);
114 
115  // what was jetB will now become the new jet
116  _bj_set_jetinfo(jetB, nn);
117 
118  } else {
119  // jet-beam recombination
120  _do_iB_recombination_step(jetA->_jets_index, diJ_min);
121  }
122 
123  // now update our nearest neighbour info and diJ table
124  // first reduce size of table
125  tail--; n--;
126  // Copy last jet contents and diJ info into position of jetA
127  *jetA = *tail;
128  diJ[jetA - head] = diJ[tail-head];
129 
130  // Initialise jetB's NN distance as well as updating it for
131  // other particles.
132  // NB: by having different loops for jetB == or != NULL we could
133  // perhaps save a few percent (usually avoid one if inside loop),
134  // but will not do it for now because on laptop fluctuations are
135  // too large to reliably measure a few percent difference...
136  for (BJ * jetI = head; jetI != tail; jetI++) {
137  // see if jetI had jetA or jetB as a NN -- if so recalculate the NN
138  if (jetI->NN == jetA || jetI->NN == jetB) {
139  _bj_set_NN_nocross(jetI, head, tail);
140  diJ[jetI-head] = _bj_diJ(jetI); // update diJ
141  }
142  // check whether new jetB is closer than jetI's current NN and
143  // if need be update things
144  if (jetB != NULL) {
145  double dist = _bj_dist(jetI,jetB);
146  if (dist < jetI->NN_dist) {
147  if (jetI != jetB) {
148  jetI->NN_dist = dist;
149  jetI->NN = jetB;
150  diJ[jetI-head] = _bj_diJ(jetI); // update diJ...
151  }
152  }
153  if (dist < jetB->NN_dist) {
154  if (jetI != jetB) {
155  jetB->NN_dist = dist;
156  jetB->NN = jetI;}
157  }
158  }
159  // if jetI's NN is the new tail then relabel it so that it becomes jetA
160  if (jetI->NN == tail) {jetI->NN = jetA;}
161  }
162 
163 
164  if (jetB != NULL) {diJ[jetB-head] = _bj_diJ(jetB);}
165 
166 
167  }
168 
169 
170  // final cleaning up;
171  delete[] diJ;
172  delete[] briefjets;
173 }
174 
175 
176 
177 // //----------------------------------------------------------------------
178 // // initialises a GenBriefJet
179 // template<> inline void ClusterSequence::_bj_set_jetinfo(
180 // GenBriefJet * const jetA, const int _jets_index) const {
181 //
182 // jetA->init(_jets[_jets_index]);
183 // jetA->_jets_index = _jets_index;
184 //
185 // }
186 //
187 //
188 // //----------------------------------------------------------------------
189 // // returns the distance between two GenBriefJets
190 // template<> double ClusterSequence::_bj_dist(
191 // const GenBriefJet * const jeta,
192 // const GenBriefJet * const jetb) const {
193 // return jeta->geom_ij(jetb);
194 // }
195 
196 FASTJET_END_NAMESPACE
197 
198 #endif // __FASTJET_CLUSTERQUENCE_N2_ICC__
199 
200