FastJet  3.3.3
NNH.hh
1 #ifndef __FASTJET_NNH_HH__
2 #define __FASTJET_NNH_HH__
3 
4 //FJSTARTHEADER
5 // $Id: NNH.hh 4420 2019-11-29 09:28:20Z soyez $
6 //
7 // Copyright (c) 2005-2019, 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 
36 FASTJET_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 ///
75 template<class BJ, class I = _NoInfo> class NNH : public NNBase<I> {
76 public:
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 
102 private:
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 //----------------------------------------------------------------------
157 template<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 //----------------------------------------------------------------------
187 template<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
209 template<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 //----------------------------------------------------------------------
230 template<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
287 template <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
310 template <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 
343 FASTJET_END_NAMESPACE // defined in fastjet/internal/base.hh
344 
345 
346 #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:80
Help solve closest pair problems with generic interparticle and beam distance (generic case) ...
Definition: NNH.hh:75
~NNH()
a destructor
Definition: NNH.hh:98
Helps solve closest pair problems with generic interparticle and particle-beam distances.
Definition: NNBase.hh:164
Class to contain pseudojets, including minimal information of use to jet-clustering routines...
Definition: PseudoJet.hh:67