Fork me on GitHub

source: git/external/fastjet/NNFJN2Plain.hh@ 0e6b14e

ImprovedOutputFile Timing
Last change on this file since 0e6b14e was b7b836a, checked in by Pavel Demin <pavel-demin@…>, 7 years ago

update FastJet library to 3.3.1 and FastJet Contrib library to 1.036

  • Property mode set to 100644
File size: 13.2 KB
Line 
1#ifndef __FASTJET_NNFJN2PLAIN_HH__
2#define __FASTJET_NNFJN2PLAIN_HH__
3
4//FJSTARTHEADER
5// $Id: NNFJN2Plain.hh 4355 2018-04-22 15:38:54Z 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 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///
97template<class BJ, class I = _NoInfo> class NNFJN2Plain : public NNBase<I> {
98public:
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
120 ~NNFJN2Plain() {
121 delete[] briefjets;
122 delete[] diJ;
123 }
124
125private:
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//----------------------------------------------------------------------
193template<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//----------------------------------------------------------------------
232template<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
254template<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//----------------------------------------------------------------------
279template<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
342template <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
365template <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
398FASTJET_END_NAMESPACE // defined in fastjet/internal/base.hh
399
400
401#endif // __FASTJET_NNFJN2PLAIN_HH__
Note: See TracBrowser for help on using the repository browser.