FastJet  3.2.1
NNFJN2Plain.hh
1 #ifndef __FASTJET_NNFJN2PLAIN_HH__
2 #define __FASTJET_NNFJN2PLAIN_HH__
3 
4 //FJSTARTHEADER
5 // $Id: NNFJN2Plain.hh 4058 2016-03-03 15:39:40Z salam $
6 //
7 // Copyright (c) 2016, 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 NNFJN2Plain
41 ///
42 /// Helps solve closest pair problems with factorised interparticle
43 /// and beam distances (ie satisfying the FastJet lemma)
44 ///
45 /// (see NNBase.hh for an introductory description)
46 ///
47 /// This variant provides an implementation based on the N2Plain
48 /// clustering strategy in FastJet. The interparticle and beam
49 /// distances should be of the form
50 ///
51 /// \code
52 /// dij = min(mom_factor(i), mom_factor(j)) * geometrical_distance(i,j)
53 /// diB = mom_factor(i) * geometrical_beam_distance(i)
54 /// \endcode
55 ///
56 /// The class is templated with a BJ (brief jet) class and can be used
57 /// with or without an extra "Information" template,
58 /// i.e. NNFJN2Plain<BJ> or NNFJN2Plain<BJ,I>
59 ///
60 /// For the NNH_N2Plain<BJ> version of the class to function, BJ must
61 /// provide four member functions
62 ///
63 /// \code
64 /// void BJ::init(const PseudoJet & jet); // initialise with a PseudoJet
65 /// double BJ::geometrical_distance(const BJ * other_bj_jet); // distance between this and other_bj_jet (geometrical part)
66 /// double BJ::geometrical_beam_distance(); // distance to the beam (geometrical part)
67 /// double BJ::momentum_factor(); // extra momentum factor
68 /// \endcode
69 ///
70 /// For the NNH_N2Plain<BJ,I> version to function, the BJ::init(...)
71 /// member must accept an extra argument
72 ///
73 /// \code
74 /// void BJ::init(const PseudoJet & jet, I * info); // initialise with a PseudoJet + info
75 /// \endcode
76 ///
77 /// NOTE: THE DISTANCE MUST BE SYMMETRIC I.E. SATISFY
78 /// \code
79 /// a.geometrical_distance(b) == b.geometrical_distance(a)
80 /// \endcode
81 ///
82 /// Note that you are strongly advised to add the following lines to
83 /// your BJ class to allow it to be used also with NNH:
84 ///
85 /// \code
86 /// /// make this BJ class compatible with the use of NNH
87 /// double BJ::distance(const BJ * other_bj_jet){
88 /// double mom1 = momentum_factor();
89 /// double mom2 = other_bj_jet->momentum_factor();
90 /// return (mom1<mom2 ? mom1 : mom2) * geometrical_distance(other_bj_jet);
91 /// }
92 /// double BJ::beam_distance(){
93 /// return momentum_factor() * geometrical_beam_distance();
94 /// }
95 /// \endcode
96 ///
97 template<class BJ, class I = _NoInfo> class NNFJN2Plain : public NNBase<I> {
98 public:
99 
100  /// constructor with an initial set of jets (which will be assigned indices
101  /// `0...jets.size()-1`)
102  NNFJN2Plain(const std::vector<PseudoJet> & jets) : NNBase<I>() {start(jets);}
103  NNFJN2Plain(const std::vector<PseudoJet> & jets, I * info) : NNBase<I>(info) {start(jets);}
104 
105  /// initialisation from a given list of particles
106  virtual void start(const std::vector<PseudoJet> & jets);
107 
108  /// returns the dij_min and indices iA, iB, for the corresponding jets.
109  /// If iB < 0 then iA recombines with the beam
110  double dij_min(int & iA, int & iB);
111 
112  /// removes the jet pointed to by index iA
113  void remove_jet(int iA);
114 
115  /// merges the jets pointed to by indices A and B and replace them with
116  /// jet, assigning it an index jet_index.
117  void merge_jets(int iA, int iB, const PseudoJet & jet, int jet_index);
118 
119  /// a destructor
121  delete[] briefjets;
122  delete[] diJ;
123  }
124 
125 private:
126  class NNBJ; // forward declaration
127 
128  // return the full distance of a particle to its NN
129  inline double compute_diJ(const NNBJ * const jet) const {
130  double mom_fact = jet->momentum_factor();
131  if (jet->NN != NULL) {
132  double other_mom_fact = jet->NN->momentum_factor();
133  if (other_mom_fact < mom_fact) {mom_fact = other_mom_fact;}
134  }
135  return jet->NN_dist * mom_fact;
136  }
137 
138  /// establish the nearest neighbour for jet, and cross check consistency
139  /// of distances for the other jets that are encountered. Assumes
140  /// jet not contained within begin...end
141  void set_NN_crosscheck(NNBJ * jet, NNBJ * begin, NNBJ * end);
142 
143  /// establish the nearest neighbour for jet; don't cross check other jets'
144  /// distances and allow jet to be contained within begin...end
145  void set_NN_nocross (NNBJ * jet, NNBJ * begin, NNBJ * end);
146 
147  /// contains the briefjets
148  NNBJ * briefjets;
149 
150  /// semaphores for the current extent of our structure
151  NNBJ * head, * tail;
152 
153  /// currently active number of jets
154  int n;
155 
156  /// where_is[i] contains a pointer to the jet with index i
157  std::vector<NNBJ *> where_is;
158 
159  /// a table containing all the (full) distances to each object's NN
160  double * diJ;
161 
162  /// a class that wraps around the BJ, supplementing it with extra information
163  /// such as pointers to neighbours, etc.
164  class NNBJ : public BJ {
165  public:
166  void init(const PseudoJet & jet, int index_in) {
167  BJ::init(jet);
168  other_init(index_in);
169  }
170  void init(const PseudoJet & jet, int index_in, I * info) {
171  BJ::init(jet, info);
172  other_init(index_in);
173  }
174  void other_init(int index_in) {
175  _index = index_in;
176  NN_dist = BJ::geometrical_beam_distance();
177  NN = NULL;
178  }
179  int index() const {return _index;}
180 
181  double NN_dist;
182  NNBJ * NN;
183 
184  private:
185  int _index;
186  };
187 
188 };
189 
190 
191 
192 //----------------------------------------------------------------------
193 template<class BJ, class I> void NNFJN2Plain<BJ,I>::start(const std::vector<PseudoJet> & jets) {
194  n = jets.size();
195  briefjets = new NNBJ[n];
196  where_is.resize(2*n);
197 
198  NNBJ * jetA = briefjets;
199 
200  // initialise the basic jet info
201  for (int i = 0; i< n; i++) {
202  // the this-> in the next line is required by standard compiler
203  // see e.g. http://stackoverflow.com/questions/10639053/name-lookups-in-c-templates
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 
220  // now create the diJ (where J is i's NN) table -- remember that
221  // we differ from standard normalisation here by a factor of R2
222  diJ = new double[n];
223  jetA = head;
224  for (int i = 0; i < n; i++) {
225  diJ[i] = compute_diJ(jetA);
226  jetA++; // have jetA follow i
227  }
228 }
229 
230 
231 //----------------------------------------------------------------------
232 template<class BJ, class I> double NNFJN2Plain<BJ,I>::dij_min(int & iA, int & iB) {
233  // find the minimum of the diJ on this round
234  double diJ_min = diJ[0];
235  int diJ_min_jet = 0;
236  for (int i = 1; i < n; i++) {
237  if (diJ[i] < diJ_min) {
238  diJ_min_jet = i;
239  diJ_min = diJ[i];
240  }
241  }
242 
243  // return information to user about recombination
244  NNBJ * jetA = & briefjets[diJ_min_jet];
245  //std::cout << jetA - briefjets << std::endl;
246  iA = jetA->index();
247  iB = jetA->NN ? jetA->NN->index() : -1;
248  return diJ_min;
249 }
250 
251 
252 //----------------------------------------------------------------------
253 // remove jetA from the list
254 template<class BJ, class I> void NNFJN2Plain<BJ,I>::remove_jet(int iA) {
255  NNBJ * jetA = where_is[iA];
256  // now update our nearest neighbour info and diJ table
257  // first reduce size of table
258  tail--; n--;
259  // Copy last jet contents and diJ info into position of jetA
260  *jetA = *tail;
261  // update the info on where the given index is stored
262  where_is[jetA->index()] = jetA;
263  diJ[jetA - head] = diJ[tail-head];
264 
265  // updating NN infos. 2 cases to watch for (see below)
266  for (NNBJ * jetI = head; jetI != tail; jetI++) {
267  // see if jetI had jetA as a NN -- if so recalculate the NN
268  if (jetI->NN == jetA) {
269  set_NN_nocross(jetI, head, tail);
270  diJ[jetI-head] = compute_diJ(jetI); // update diJ
271  }
272  // if jetI's NN is the new tail then relabel it so that it becomes jetA
273  if (jetI->NN == tail) {jetI->NN = jetA;}
274  }
275 }
276 
277 
278 //----------------------------------------------------------------------
279 template<class BJ, class I> void NNFJN2Plain<BJ,I>::merge_jets(int iA, int iB,
280  const PseudoJet & jet, int index) {
281 
282  NNBJ * jetA = where_is[iA];
283  NNBJ * jetB = where_is[iB];
284 
285  // If necessary relabel A & B to ensure jetB < jetA, that way if
286  // the larger of them == newtail then that ends up being jetA and
287  // the new jet that is added as jetB is inserted in a position that
288  // has a future!
289  if (jetA < jetB) std::swap(jetA,jetB);
290 
291  // initialise jetB based on the new jet
292  this->init_jet(jetB, jet, index);
293  // and record its position (making sure we have the space)
294  if (index >= int(where_is.size())) where_is.resize(2*index);
295  where_is[jetB->index()] = jetB;
296 
297  // now update our nearest neighbour info
298  // first reduce size of table
299  tail--; n--;
300  // Copy last jet contents into position of jetA
301  *jetA = *tail;
302  // update the info on where the tail's index is stored
303  where_is[jetA->index()] = jetA;
304  diJ[jetA - head] = diJ[tail-head];
305 
306  // initialise jetB NN's distance and update NN infos
307  for (NNBJ * jetI = head; jetI != tail; jetI++) {
308  // see if jetI had jetA or jetB as a NN -- if so recalculate the NN
309  if (jetI->NN == jetA || jetI->NN == jetB) {
310  set_NN_nocross(jetI, head, tail);
311  diJ[jetI-head] = compute_diJ(jetI); // update diJ
312  }
313 
314  // check whether new jetB is closer than jetI's current NN and
315  // if need be update things
316  double dist = jetI->geometrical_distance(jetB);
317  if (dist < jetI->NN_dist) { // we need to update I
318  if (jetI != jetB) {
319  jetI->NN_dist = dist;
320  jetI->NN = jetB;
321  diJ[jetI-head] = compute_diJ(jetI); // update diJ...
322  }
323  }
324  if (dist < jetB->NN_dist) { // we need to update B
325  if (jetI != jetB) {
326  jetB->NN_dist = dist;
327  jetB->NN = jetI;
328  }
329  }
330 
331  // if jetI's NN is the new tail then relabel it so that it becomes jetA
332  if (jetI->NN == tail) {jetI->NN = jetA;}
333  }
334 
335  // update info for jetB
336  diJ[jetB-head] = compute_diJ(jetB);
337 }
338 
339 
340 //----------------------------------------------------------------------
341 // this function assumes that jet is not contained within begin...end
342 template <class BJ, class I> void NNFJN2Plain<BJ,I>::set_NN_crosscheck(NNBJ * jet,
343  NNBJ * begin, NNBJ * end) {
344  double NN_dist = jet->geometrical_beam_distance();
345  NNBJ * NN = NULL;
346  for (NNBJ * jetB = begin; jetB != end; jetB++) {
347  double dist = jet->geometrical_distance(jetB);
348  if (dist < NN_dist) {
349  NN_dist = dist;
350  NN = jetB;
351  }
352  if (dist < jetB->NN_dist) {
353  jetB->NN_dist = dist;
354  jetB->NN = jet;
355  }
356  }
357  jet->NN = NN;
358  jet->NN_dist = NN_dist;
359 }
360 
361 
362 //----------------------------------------------------------------------
363 // set the NN for jet without checking whether in the process you might
364 // have discovered a new nearest neighbour for another jet
365 template <class BJ, class I> void NNFJN2Plain<BJ,I>::set_NN_nocross(
366  NNBJ * jet, NNBJ * begin, NNBJ * end) {
367  double NN_dist = jet->geometrical_beam_distance();
368  NNBJ * NN = NULL;
369  // if (head < jet) {
370  // for (NNBJ * jetB = head; jetB != jet; jetB++) {
371  if (begin < jet) {
372  for (NNBJ * jetB = begin; jetB != jet; jetB++) {
373  double dist = jet->geometrical_distance(jetB);
374  if (dist < NN_dist) {
375  NN_dist = dist;
376  NN = jetB;
377  }
378  }
379  }
380  // if (tail > jet) {
381  // for (NNBJ * jetB = jet+1; jetB != tail; jetB++) {
382  if (end > jet) {
383  for (NNBJ * jetB = jet+1; jetB != end; jetB++) {
384  double dist = jet->geometrical_distance(jetB);
385  if (dist < NN_dist) {
386  NN_dist = dist;
387  NN = jetB;
388  }
389  }
390  }
391  jet->NN = NN;
392  jet->NN_dist = NN_dist;
393 }
394 
395 
396 
397 
398 FASTJET_END_NAMESPACE // defined in fastjet/internal/base.hh
399 
400 
401 #endif // __FASTJET_NNFJN2PLAIN_HH__
~NNFJN2Plain()
a destructor
Definition: NNFJN2Plain.hh:120
Helps solve closest pair problems with generic interparticle and particle-beam distances.
Definition: NNBase.hh:164
NNFJN2Plain(const std::vector< PseudoJet > &jets)
constructor with an initial set of jets (which will be assigned indices 0...jets.size()-1) ...
Definition: NNFJN2Plain.hh:102
Class to contain pseudojets, including minimal information of use to jet-clustering routines...
Definition: PseudoJet.hh:67
Helps solve closest pair problems with factorised interparticle and beam distances (ie satisfying the...
Definition: NNFJN2Plain.hh:97