FastJet  3.1.3
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Modules Pages
NNH.hh
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 
37 FASTJET_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
42 class _NoInfo {};
43 
44 /// @ingroup advanced_usage
45 /// \class NNHInfo
46 /// template that will help initialise a BJ with a PseudoJet and extra information
47 template<class I> class NNHInfo {
48 public:
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);}
52 private:
53  I * _info;
54 };
55 
56 /// @ingroup advanced_usage
57 /// Specialisation of NNHInfo for cases where there is no extra info
58 template<> class NNHInfo<_NoInfo> {
59 public:
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 
113 template<class BJ, class I = _NoInfo> class NNH : public NNHInfo<I> {
114 public:
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 
139 private:
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 //----------------------------------------------------------------------
194 template<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 //----------------------------------------------------------------------
224 template<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
246 template<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 //----------------------------------------------------------------------
267 template<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
324 template <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
347 template <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 
380 FASTJET_END_NAMESPACE // defined in fastjet/internal/base.hh
381 
382 
383 #endif // __FASTJET_NNH_HH__
NNH(const std::vector< PseudoJet > &jets)
constructor with an initial set of jets (which will be assigned indices 0 ...
Definition: NNH.hh:118
Help solve closest pair problems with generic interparticle and beam distance.
Definition: NNH.hh:113
~NNH()
a destructor
Definition: NNH.hh:135
template that will help initialise a BJ with a PseudoJet and extra information
Definition: NNH.hh:47
dummy class, used as a default template argument
Definition: NNH.hh:42
Class to contain pseudojets, including minimal information of use to jet-clustering routines...
Definition: PseudoJet.hh:67