Fork me on GitHub

source: git/external/fastjet/NNH.hh@ 1d5395b

Timing
Last change on this file since 1d5395b was b7b836a, checked in by Pavel Demin <pavel-demin@…>, 6 years ago

update FastJet library to 3.3.1 and FastJet Contrib library to 1.036

  • Property mode set to 100644
File size: 11.1 KB
Line 
1#ifndef __FASTJET_NNH_HH__
2#define __FASTJET_NNH_HH__
3
4//FJSTARTHEADER
5// $Id: NNH.hh 4354 2018-04-22 07:12:37Z salam $
6//
7// Copyright (c) 2005-2018, 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/NNBase.hh>
35
36FASTJET_BEGIN_NAMESPACE // defined in fastjet/internal/base.hh
37
38//----------------------------------------------------------------------
39/// @ingroup advanced_usage
40/// \class NNH
41/// Help solve closest pair problems with generic interparticle and
42/// beam distance (generic case)
43///
44/// (see NNBase.hh for an introductory description)
45///
46/// This variant provides an implementation for any distance measure.
47/// It is templated with a BJ (brief jet) classand can be used with or
48/// without an extra "Information" template, i.e. NNH<BJ> or NNH<BJ,I>
49///
50/// For the NNH<BJ> version of the class to function, BJ must provide
51/// three member functions
52///
53/// - void BJ::init(const PseudoJet & jet); // initialise with a PseudoJet
54/// - double BJ::distance(const BJ * other_bj_jet); // distance between this and other_bj_jet
55/// - double BJ::beam_distance() ; // distance to the beam
56///
57/// For the NNH<BJ,I> version to function, the BJ::init(...) member
58/// must accept an extra argument
59///
60/// - void BJ::init(const PseudoJet & jet, I * info); // initialise with a PseudoJet + info
61///
62/// NOTE: THE DISTANCE MUST BE SYMMETRIC I.E. SATISFY
63/// a.distance(b) == b.distance(a)
64///
65/// For an example of how the NNH<BJ> class is used, see the Jade (and
66/// EECambridge) plugins
67///
68/// NB: the NNH algorithm is expected N^2, but has a worst case of
69/// N^3. Many QCD problems tend to place one closer to the N^3 end of
70/// the spectrum than one would like. There is scope for further
71/// progress (cf Eppstein, Cardinal & Eppstein), nevertheless the
72/// current class is already significantly faster than standard N^3
73/// implementations.
74///
75template<class BJ, class I = _NoInfo> class NNH : public NNBase<I> {
76public:
77
78 /// constructor with an initial set of jets (which will be assigned indices
79 /// 0 ... jets.size()-1
80 NNH(const std::vector<PseudoJet> & jets) : NNBase<I>() {start(jets);}
81 NNH(const std::vector<PseudoJet> & jets, I * info) : NNBase<I>(info) {start(jets);}
82
83 // initialisation from a given list of particles
84 void start(const std::vector<PseudoJet> & jets);
85
86 /// return the dij_min and indices iA, iB, for the corresponding jets.
87 /// If iB < 0 then iA recombines with the beam
88 double dij_min(int & iA, int & iB);
89
90 /// remove the jet pointed to by index iA
91 void remove_jet(int iA);
92
93 /// merge the jets pointed to by indices A and B and replace them with
94 /// jet, assigning it an index jet_index.
95 void merge_jets(int iA, int iB, const PseudoJet & jet, int jet_index);
96
97 /// a destructor
98 ~NNH() {
99 delete[] briefjets;
100 }
101
102private:
103 class NNBJ; // forward declaration
104
105 /// establish the nearest neighbour for jet, and cross check consistency
106 /// of distances for the other jets that are encountered. Assumes
107 /// jet not contained within begin...end
108 void set_NN_crosscheck(NNBJ * jet, NNBJ * begin, NNBJ * end);
109
110 /// establish the nearest neighbour for jet; don't cross check other jets'
111 /// distances and allow jet to be contained within begin...end
112 void set_NN_nocross (NNBJ * jet, NNBJ * begin, NNBJ * end);
113
114 /// contains the briefjets
115 NNBJ * briefjets;
116
117 /// semaphores for the current extent of our structure
118 NNBJ * head, * tail;
119
120 /// currently active number of jets
121 int n;
122
123 /// where_is[i] contains a pointer to the jet with index i
124 std::vector<NNBJ *> where_is;
125
126 /// a class that wraps around the BJ, supplementing it with extra information
127 /// such as pointers to neighbours, etc.
128 class NNBJ : public BJ {
129 public:
130 void init(const PseudoJet & jet, int index_in) {
131 BJ::init(jet);
132 other_init(index_in);
133 }
134 void init(const PseudoJet & jet, int index_in, I * info) {
135 BJ::init(jet, info);
136 other_init(index_in);
137 }
138 void other_init(int index_in) {
139 _index = index_in;
140 NN_dist = BJ::beam_distance();
141 NN = NULL;
142 }
143 int index() const {return _index;}
144
145 double NN_dist;
146 NNBJ * NN;
147
148 private:
149 int _index;
150 };
151
152};
153
154
155
156//----------------------------------------------------------------------
157template<class BJ, class I> void NNH<BJ,I>::start(const std::vector<PseudoJet> & jets) {
158 n = jets.size();
159 briefjets = new NNBJ[n];
160 where_is.resize(2*n);
161
162 NNBJ * jetA = briefjets;
163
164 // initialise the basic jet info
165 for (int i = 0; i< n; i++) {
166 //jetA->init(jets[i], i);
167 this->init_jet(jetA, jets[i], i);
168 where_is[i] = jetA;
169 jetA++; // move on to next entry of briefjets
170 }
171 tail = jetA; // a semaphore for the end of briefjets
172 head = briefjets; // a nicer way of naming start
173
174 // now initialise the NN distances: jetA will run from 1..n-1; and
175 // jetB from 0..jetA-1
176 for (jetA = head + 1; jetA != tail; jetA++) {
177 // set NN info for jetA based on jets running from head..jetA-1,
178 // checking in the process whether jetA itself is an undiscovered
179 // NN of one of those jets.
180 set_NN_crosscheck(jetA, head, jetA);
181 }
182 //std::cout << "OOOO " << briefjets[1].NN_dist << " " << briefjets[1].NN - briefjets << std::endl;
183}
184
185
186//----------------------------------------------------------------------
187template<class BJ, class I> double NNH<BJ,I>::dij_min(int & iA, int & iB) {
188 // find the minimum of the diJ on this round
189 double diJ_min = briefjets[0].NN_dist;
190 int diJ_min_jet = 0;
191 for (int i = 1; i < n; i++) {
192 if (briefjets[i].NN_dist < diJ_min) {
193 diJ_min_jet = i;
194 diJ_min = briefjets[i].NN_dist;
195 }
196 }
197
198 // return information to user about recombination
199 NNBJ * jetA = & briefjets[diJ_min_jet];
200 //std::cout << jetA - briefjets << std::endl;
201 iA = jetA->index();
202 iB = jetA->NN ? jetA->NN->index() : -1;
203 return diJ_min;
204}
205
206
207//----------------------------------------------------------------------
208// remove jetA from the list
209template<class BJ, class I> void NNH<BJ,I>::remove_jet(int iA) {
210 NNBJ * jetA = where_is[iA];
211 // now update our nearest neighbour info and diJ table
212 // first reduce size of table
213 tail--; n--;
214 // Copy last jet contents and diJ info into position of jetA
215 *jetA = *tail;
216 // update the info on where the given index is stored
217 where_is[jetA->index()] = jetA;
218
219 for (NNBJ * jetI = head; jetI != tail; jetI++) {
220 // see if jetI had jetA or jetB as a NN -- if so recalculate the NN
221 if (jetI->NN == jetA) set_NN_nocross(jetI, head, tail);
222
223 // if jetI's NN is the new tail then relabel it so that it becomes jetA
224 if (jetI->NN == tail) {jetI->NN = jetA;}
225 }
226}
227
228
229//----------------------------------------------------------------------
230template<class BJ, class I> void NNH<BJ,I>::merge_jets(int iA, int iB,
231 const PseudoJet & jet, int index) {
232
233 NNBJ * jetA = where_is[iA];
234 NNBJ * jetB = where_is[iB];
235
236 // If necessary relabel A & B to ensure jetB < jetA, that way if
237 // the larger of them == newtail then that ends up being jetA and
238 // the new jet that is added as jetB is inserted in a position that
239 // has a future!
240 if (jetA < jetB) std::swap(jetA,jetB);
241
242 // initialise jetB based on the new jet
243 //jetB->init(jet, index);
244 this->init_jet(jetB, jet, index);
245 // and record its position (making sure we have the space)
246 if (index >= int(where_is.size())) where_is.resize(2*index);
247 where_is[jetB->index()] = jetB;
248
249 // now update our nearest neighbour info
250 // first reduce size of table
251 tail--; n--;
252 // Copy last jet contents into position of jetA
253 *jetA = *tail;
254 // update the info on where the tail's index is stored
255 where_is[jetA->index()] = jetA;
256
257 for (NNBJ * jetI = head; jetI != tail; jetI++) {
258 // see if jetI had jetA or jetB as a NN -- if so recalculate the NN
259 if (jetI->NN == jetA || jetI->NN == jetB) {
260 set_NN_nocross(jetI, head, tail);
261 }
262
263 // check whether new jetB is closer than jetI's current NN and
264 // if need be update things
265 double dist = jetI->distance(jetB);
266 if (dist < jetI->NN_dist) {
267 if (jetI != jetB) {
268 jetI->NN_dist = dist;
269 jetI->NN = jetB;
270 }
271 }
272 if (dist < jetB->NN_dist) {
273 if (jetI != jetB) {
274 jetB->NN_dist = dist;
275 jetB->NN = jetI;
276 }
277 }
278
279 // if jetI's NN is the new tail then relabel it so that it becomes jetA
280 if (jetI->NN == tail) {jetI->NN = jetA;}
281 }
282}
283
284
285//----------------------------------------------------------------------
286// this function assumes that jet is not contained within begin...end
287template <class BJ, class I> void NNH<BJ,I>::set_NN_crosscheck(NNBJ * jet,
288 NNBJ * begin, NNBJ * end) {
289 double NN_dist = jet->beam_distance();
290 NNBJ * NN = NULL;
291 for (NNBJ * jetB = begin; jetB != end; jetB++) {
292 double dist = jet->distance(jetB);
293 if (dist < NN_dist) {
294 NN_dist = dist;
295 NN = jetB;
296 }
297 if (dist < jetB->NN_dist) {
298 jetB->NN_dist = dist;
299 jetB->NN = jet;
300 }
301 }
302 jet->NN = NN;
303 jet->NN_dist = NN_dist;
304}
305
306
307//----------------------------------------------------------------------
308// set the NN for jet without checking whether in the process you might
309// have discovered a new nearest neighbour for another jet
310template <class BJ, class I> void NNH<BJ,I>::set_NN_nocross(
311 NNBJ * jet, NNBJ * begin, NNBJ * end) {
312 double NN_dist = jet->beam_distance();
313 NNBJ * NN = NULL;
314 // if (head < jet) {
315 // for (NNBJ * jetB = head; jetB != jet; jetB++) {
316 if (begin < jet) {
317 for (NNBJ * jetB = begin; jetB != jet; jetB++) {
318 double dist = jet->distance(jetB);
319 if (dist < NN_dist) {
320 NN_dist = dist;
321 NN = jetB;
322 }
323 }
324 }
325 // if (tail > jet) {
326 // for (NNBJ * jetB = jet+1; jetB != tail; jetB++) {
327 if (end > jet) {
328 for (NNBJ * jetB = jet+1; jetB != end; jetB++) {
329 double dist = jet->distance (jetB);
330 if (dist < NN_dist) {
331 NN_dist = dist;
332 NN = jetB;
333 }
334 }
335 }
336 jet->NN = NN;
337 jet->NN_dist = NN_dist;
338}
339
340
341
342
343FASTJET_END_NAMESPACE // defined in fastjet/internal/base.hh
344
345
346#endif // __FASTJET_NNH_HH__
Note: See TracBrowser for help on using the repository browser.