Fork me on GitHub

source: git/external/fastjet/NNH.hh@ d091310

ImprovedOutputFile Timing dual_readout llp
Last change on this file since d091310 was 35cdc46, checked in by Pavel Demin <demin@…>, 10 years ago

upgrade FastJet to version 3.1.0-beta.1, upgrade Nsubjettiness to version 2.1.0, add SoftKiller version 1.0.0

  • Property mode set to 100644
File size: 12.3 KB
Line 
1#ifndef __FASTJET_NNH_HH__
2#define __FASTJET_NNH_HH__
3
4//FJSTARTHEADER
5// $Id: NNH.hh 3433 2014-07-23 08:17:03Z salam $
6//
7// Copyright (c) 2005-2014, Matteo Cacciari, Gavin P. Salam and Gregory Soyez
8//
9//----------------------------------------------------------------------
10// This file is part of FastJet.
11//
12// FastJet is free software; you can redistribute it and/or modify
13// it under the terms of the GNU General Public License as published by
14// the Free Software Foundation; either version 2 of the License, or
15// (at your option) any later version.
16//
17// The algorithms that underlie FastJet have required considerable
18// development. They are described in the original FastJet paper,
19// hep-ph/0512210 and in the manual, arXiv:1111.6097. If you use
20// FastJet as part of work towards a scientific publication, please
21// quote the version you use and include a citation to the manual and
22// optionally also to hep-ph/0512210.
23//
24// FastJet is distributed in the hope that it will be useful,
25// but WITHOUT ANY WARRANTY; without even the implied warranty of
26// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
27// GNU General Public License for more details.
28//
29// You should have received a copy of the GNU General Public License
30// along with FastJet. If not, see <http://www.gnu.org/licenses/>.
31//----------------------------------------------------------------------
32//FJENDHEADER
33
34#include<fastjet/ClusterSequence.hh>
35
36
37FASTJET_BEGIN_NAMESPACE // defined in fastjet/internal/base.hh
38
39/// @ingroup advanced_usage
40/// \class _NoInfo
41/// dummy class, used as a default template argument
42class _NoInfo {};
43
44/// @ingroup advanced_usage
45/// \class NNHInfo
46/// template that will help initialise a BJ with a PseudoJet and extra information
47template<class I> class NNHInfo {
48public:
49 NNHInfo() : _info(NULL) {}
50 NNHInfo(I * info) : _info(info) {}
51 template<class NNBJ> void init_jet(NNBJ * briefjet, const fastjet::PseudoJet & jet, int index) { briefjet->init(jet, index, _info);}
52private:
53 I * _info;
54};
55
56/// @ingroup advanced_usage
57/// Specialisation of NNHInfo for cases where there is no extra info
58template<> class NNHInfo<_NoInfo> {
59public:
60 NNHInfo() {}
61 NNHInfo(_NoInfo * ) {}
62 template<class NNBJ> void init_jet(NNBJ * briefjet, const fastjet::PseudoJet & jet, int index) { briefjet->init(jet, index);}
63};
64
65
66//----------------------------------------------------------------------
67/// @ingroup advanced_usage
68/// \class NNH
69/// Help solve closest pair problems with generic interparticle and
70/// beam distance.
71///
72/// Class to help solve closest pair problems with generic interparticle
73/// distances and a beam distance, using Anderberg's Nearest Neighbour
74/// Heuristic.
75///
76/// It is templated with a BJ (brief jet) class --- BJ should
77/// basically cache the minimal amount of information that is needed
78/// to efficiently calculate interparticle distances and particle-beam
79/// distances.
80///
81/// This class can be used with or without an extra "Information" template,
82/// i.e. NNB<BJ> or NNH<BJ,I>
83///
84/// For the NNH<BJ> version of the class to function, BJ must provide
85/// three member functions
86///
87/// - void BJ::init(const PseudoJet & jet); // initialise with a PseudoJet
88/// - double BJ::distance(const BJ * other_bj_jet); // distance between this and other_bj_jet
89/// - double BJ::beam_distance() ; // distance to the beam
90///
91/// For the NNH<BJ,I> version to function, the BJ::init(...) member
92/// must accept an extra argument
93///
94/// - void BJ::init(const PseudoJet & jet, I * info); // initialise with a PseudoJet + info
95///
96/// where info might be a pointer to a class that contains, e.g., information
97/// about R, or other parameters of the jet algorithm
98///
99/// For an example of how the NNH<BJ> class is used, see the Jade (and
100/// EECambridge) plugins
101///
102/// NB: the NNH algorithm is expected N^2, but has a worst case of
103/// N^3. Many QCD problems tend to place one closer to the N^3 end of
104/// the spectrum than one would like. There is scope for further
105/// progress (cf Eppstein, Cardinal & Eppstein), nevertheless the
106/// current class is already significantly faster than standard N^3
107/// implementations.
108///
109///
110/// Implementation note: this class derives from NNHInfo, which deals
111/// with storing any global information that is needed during the clustering
112
113template<class BJ, class I = _NoInfo> class NNH : public NNHInfo<I> {
114public:
115
116 /// constructor with an initial set of jets (which will be assigned indices
117 /// 0 ... jets.size()-1
118 NNH(const std::vector<PseudoJet> & jets) {start(jets);}
119 NNH(const std::vector<PseudoJet> & jets, I * info) : NNHInfo<I>(info) {start(jets);}
120
121 void start(const std::vector<PseudoJet> & jets);
122
123 /// return the dij_min and indices iA, iB, for the corresponding jets.
124 /// If iB < 0 then iA recombines with the beam
125 double dij_min(int & iA, int & iB);
126
127 /// remove the jet pointed to by index iA
128 void remove_jet(int iA);
129
130 /// merge the jets pointed to by indices A and B and replace them with
131 /// jet, assigning it an index jet_index.
132 void merge_jets(int iA, int iB, const PseudoJet & jet, int jet_index);
133
134 /// a destructor
135 ~NNH() {
136 delete[] briefjets;
137 }
138
139private:
140 class NNBJ; // forward declaration
141
142 /// establish the nearest neighbour for jet, and cross check consistency
143 /// of distances for the other jets that are encountered. Assumes
144 /// jet not contained within begin...end
145 void set_NN_crosscheck(NNBJ * jet, NNBJ * begin, NNBJ * end);
146
147 /// establish the nearest neighbour for jet; don't cross check other jets'
148 /// distances and allow jet to be contained within begin...end
149 void set_NN_nocross (NNBJ * jet, NNBJ * begin, NNBJ * end);
150
151 /// contains the briefjets
152 NNBJ * briefjets;
153
154 /// semaphores for the current extent of our structure
155 NNBJ * head, * tail;
156
157 /// currently active number of jets
158 int n;
159
160 /// where_is[i] contains a pointer to the jet with index i
161 std::vector<NNBJ *> where_is;
162
163 /// a class that wraps around the BJ, supplementing it with extra information
164 /// such as pointers to neighbours, etc.
165 class NNBJ : public BJ {
166 public:
167 void init(const PseudoJet & jet, int index_in) {
168 BJ::init(jet);
169 other_init(index_in);
170 }
171 void init(const PseudoJet & jet, int index_in, I * info) {
172 BJ::init(jet, info);
173 other_init(index_in);
174 }
175 void other_init(int index_in) {
176 _index = index_in;
177 NN_dist = BJ::beam_distance();
178 NN = NULL;
179 }
180 int index() const {return _index;}
181
182 double NN_dist;
183 NNBJ * NN;
184
185 private:
186 int _index;
187 };
188
189};
190
191
192
193//----------------------------------------------------------------------
194template<class BJ, class I> void NNH<BJ,I>::start(const std::vector<PseudoJet> & jets) {
195 n = jets.size();
196 briefjets = new NNBJ[n];
197 where_is.resize(2*n);
198
199 NNBJ * jetA = briefjets;
200
201 // initialise the basic jet info
202 for (int i = 0; i< n; i++) {
203 //jetA->init(jets[i], i);
204 this->init_jet(jetA, jets[i], i);
205 where_is[i] = jetA;
206 jetA++; // move on to next entry of briefjets
207 }
208 tail = jetA; // a semaphore for the end of briefjets
209 head = briefjets; // a nicer way of naming start
210
211 // now initialise the NN distances: jetA will run from 1..n-1; and
212 // jetB from 0..jetA-1
213 for (jetA = head + 1; jetA != tail; jetA++) {
214 // set NN info for jetA based on jets running from head..jetA-1,
215 // checking in the process whether jetA itself is an undiscovered
216 // NN of one of those jets.
217 set_NN_crosscheck(jetA, head, jetA);
218 }
219 //std::cout << "OOOO " << briefjets[1].NN_dist << " " << briefjets[1].NN - briefjets << std::endl;
220}
221
222
223//----------------------------------------------------------------------
224template<class BJ, class I> double NNH<BJ,I>::dij_min(int & iA, int & iB) {
225 // find the minimum of the diJ on this round
226 double diJ_min = briefjets[0].NN_dist;
227 int diJ_min_jet = 0;
228 for (int i = 1; i < n; i++) {
229 if (briefjets[i].NN_dist < diJ_min) {
230 diJ_min_jet = i;
231 diJ_min = briefjets[i].NN_dist;
232 }
233 }
234
235 // return information to user about recombination
236 NNBJ * jetA = & briefjets[diJ_min_jet];
237 //std::cout << jetA - briefjets << std::endl;
238 iA = jetA->index();
239 iB = jetA->NN ? jetA->NN->index() : -1;
240 return diJ_min;
241}
242
243
244//----------------------------------------------------------------------
245// remove jetA from the list
246template<class BJ, class I> void NNH<BJ,I>::remove_jet(int iA) {
247 NNBJ * jetA = where_is[iA];
248 // now update our nearest neighbour info and diJ table
249 // first reduce size of table
250 tail--; n--;
251 // Copy last jet contents and diJ info into position of jetA
252 *jetA = *tail;
253 // update the info on where the given index is stored
254 where_is[jetA->index()] = jetA;
255
256 for (NNBJ * jetI = head; jetI != tail; jetI++) {
257 // see if jetI had jetA or jetB as a NN -- if so recalculate the NN
258 if (jetI->NN == jetA) set_NN_nocross(jetI, head, tail);
259
260 // if jetI's NN is the new tail then relabel it so that it becomes jetA
261 if (jetI->NN == tail) {jetI->NN = jetA;}
262 }
263}
264
265
266//----------------------------------------------------------------------
267template<class BJ, class I> void NNH<BJ,I>::merge_jets(int iA, int iB,
268 const PseudoJet & jet, int index) {
269
270 NNBJ * jetA = where_is[iA];
271 NNBJ * jetB = where_is[iB];
272
273 // If necessary relabel A & B to ensure jetB < jetA, that way if
274 // the larger of them == newtail then that ends up being jetA and
275 // the new jet that is added as jetB is inserted in a position that
276 // has a future!
277 if (jetA < jetB) std::swap(jetA,jetB);
278
279 // initialise jetB based on the new jet
280 //jetB->init(jet, index);
281 this->init_jet(jetB, jet, index);
282 // and record its position (making sure we have the space)
283 if (index >= int(where_is.size())) where_is.resize(2*index);
284 where_is[jetB->index()] = jetB;
285
286 // now update our nearest neighbour info
287 // first reduce size of table
288 tail--; n--;
289 // Copy last jet contents into position of jetA
290 *jetA = *tail;
291 // update the info on where the tail's index is stored
292 where_is[jetA->index()] = jetA;
293
294 for (NNBJ * jetI = head; jetI != tail; jetI++) {
295 // see if jetI had jetA or jetB as a NN -- if so recalculate the NN
296 if (jetI->NN == jetA || jetI->NN == jetB) {
297 set_NN_nocross(jetI, head, tail);
298 }
299
300 // check whether new jetB is closer than jetI's current NN and
301 // if need be update things
302 double dist = jetI->distance(jetB);
303 if (dist < jetI->NN_dist) {
304 if (jetI != jetB) {
305 jetI->NN_dist = dist;
306 jetI->NN = jetB;
307 }
308 }
309 if (dist < jetB->NN_dist) {
310 if (jetI != jetB) {
311 jetB->NN_dist = dist;
312 jetB->NN = jetI;
313 }
314 }
315
316 // if jetI's NN is the new tail then relabel it so that it becomes jetA
317 if (jetI->NN == tail) {jetI->NN = jetA;}
318 }
319}
320
321
322//----------------------------------------------------------------------
323// this function assumes that jet is not contained within begin...end
324template <class BJ, class I> void NNH<BJ,I>::set_NN_crosscheck(NNBJ * jet,
325 NNBJ * begin, NNBJ * end) {
326 double NN_dist = jet->beam_distance();
327 NNBJ * NN = NULL;
328 for (NNBJ * jetB = begin; jetB != end; jetB++) {
329 double dist = jet->distance(jetB);
330 if (dist < NN_dist) {
331 NN_dist = dist;
332 NN = jetB;
333 }
334 if (dist < jetB->NN_dist) {
335 jetB->NN_dist = dist;
336 jetB->NN = jet;
337 }
338 }
339 jet->NN = NN;
340 jet->NN_dist = NN_dist;
341}
342
343
344//----------------------------------------------------------------------
345// set the NN for jet without checking whether in the process you might
346// have discovered a new nearest neighbour for another jet
347template <class BJ, class I> void NNH<BJ,I>::set_NN_nocross(
348 NNBJ * jet, NNBJ * begin, NNBJ * end) {
349 double NN_dist = jet->beam_distance();
350 NNBJ * NN = NULL;
351 // if (head < jet) {
352 // for (NNBJ * jetB = head; jetB != jet; jetB++) {
353 if (begin < jet) {
354 for (NNBJ * jetB = begin; jetB != jet; jetB++) {
355 double dist = jet->distance(jetB);
356 if (dist < NN_dist) {
357 NN_dist = dist;
358 NN = jetB;
359 }
360 }
361 }
362 // if (tail > jet) {
363 // for (NNBJ * jetB = jet+1; jetB != tail; jetB++) {
364 if (end > jet) {
365 for (NNBJ * jetB = jet+1; jetB != end; jetB++) {
366 double dist = jet->distance (jetB);
367 if (dist < NN_dist) {
368 NN_dist = dist;
369 NN = jetB;
370 }
371 }
372 }
373 jet->NN = NN;
374 jet->NN_dist = NN_dist;
375}
376
377
378
379
380FASTJET_END_NAMESPACE // defined in fastjet/internal/base.hh
381
382
383#endif // __FASTJET_NNH_HH__
Note: See TracBrowser for help on using the repository browser.