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