Analysis Software
Documentation for sPHENIX simulation software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
fjcore.cc
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file fjcore.cc
1 // fjcore -- extracted from FastJet v3.2.1 (http://fastjet.fr)
2 //
3 // fjcore constitutes a digest of the main FastJet functionality.
4 // The files fjcore.hh and fjcore.cc are meant to provide easy access to these
5 // core functions, in the form of single files and without the need of a full
6 // FastJet installation:
7 //
8 // g++ main.cc fjcore.cc
9 //
10 // with main.cc including fjcore.hh.
11 //
12 // A fortran interface, fjcorefortran.cc, is also provided. See the example
13 // and the Makefile for instructions.
14 //
15 // The results are expected to be identical to those obtained by linking to
16 // the full FastJet distribution.
17 //
18 // NOTE THAT, IN ORDER TO MAKE IT POSSIBLE FOR FJCORE AND THE FULL FASTJET
19 // TO COEXIST, THE FORMER USES THE "fjcore" NAMESPACE INSTEAD OF "fastjet".
20 //
21 // In particular, fjcore provides:
22 //
23 // - access to all native pp and ee algorithms, kt, anti-kt, C/A.
24 // For C/A, the NlnN method is available, while anti-kt and kt
25 // are limited to the N^2 one (still the fastest for N < 100k particles)
26 // - access to selectors, for implementing cuts and selections
27 // - access to all functionalities related to pseudojets (e.g. a jet's
28 // structure or user-defined information)
29 //
30 // Instead, it does NOT provide:
31 //
32 // - jet areas functionality
33 // - background estimation
34 // - access to other algorithms via plugins
35 // - interface to CGAL
36 // - fastjet tools, e.g. filters, taggers
37 //
38 // If these functionalities are needed, the full FastJet installation must be
39 // used. The code will be fully compatible, with the sole replacement of the
40 // header files and of the fjcore namespace with the fastjet one.
41 //
42 // fjcore.hh and fjcore.cc are not meant to be human-readable.
43 // For documentation, see the full FastJet manual and doxygen at http://fastjet.fr
44 //
45 // Like FastJet, fjcore is released under the terms of the GNU General Public
46 // License version 2 (GPLv2). If you use this code as part of work towards a
47 // scientific publication, whether directly or contained within another program
48 // (e.g. Delphes, MadGraph, SpartyJet, Rivet, LHC collaboration software frameworks,
49 // etc.), you should include a citation to
50 //
51 // EPJC72(2012)1896 [arXiv:1111.6097] (FastJet User Manual)
52 // and, optionally, Phys.Lett.B641 (2006) 57 [arXiv:hep-ph/0512210]
53 //
54 // Copyright (c) 2005-2016, Matteo Cacciari, Gavin P. Salam and Gregory Soyez
55 //
56 //----------------------------------------------------------------------
57 // This file is part of FastJet.
58 //
59 // FastJet is free software; you can redistribute it and/or modify
60 // it under the terms of the GNU General Public License as published by
61 // the Free Software Foundation; either version 2 of the License, or
62 // (at your option) any later version.
63 //
64 // The algorithms that underlie FastJet have required considerable
65 // development and are described in hep-ph/0512210. If you use
66 // FastJet as part of work towards a scientific publication, please
67 // include a citation to the FastJet paper.
68 //
69 // FastJet is distributed in the hope that it will be useful,
70 // but WITHOUT ANY WARRANTY; without even the implied warranty of
71 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
72 // GNU General Public License for more details.
73 //
74 // You should have received a copy of the GNU General Public License
75 // along with FastJet. If not, see <http://www.gnu.org/licenses/>.
76 //----------------------------------------------------------------------
77 //
78 #include "fjcore.hh"
79 #ifndef __FJCORE_VERSION_HH__
80 #define __FJCORE_VERSION_HH__
81 #include<string>
85 #endif // __FJCORE_VERSION_HH__
86 #ifndef __FJCORE_CLUSTERQUENCE_N2_ICC__
87 #define __FJCORE_CLUSTERQUENCE_N2_ICC__
88 FJCORE_BEGIN_NAMESPACE // defined in fastjet/internal/base.hh
89 template<class BJ> void ClusterSequence::_simple_N2_cluster() {
90  int n = _jets.size();
91  BJ * briefjets = new BJ[n];
92  BJ * jetA = briefjets, * jetB;
93  for (int i = 0; i< n; i++) {
94  _bj_set_jetinfo(jetA, i);
95  jetA++; // move on to next entry of briefjets
96  }
97  BJ * tail = jetA; // a semaphore for the end of briefjets
98  BJ * head = briefjets; // a nicer way of naming start
99  for (jetA = head + 1; jetA != tail; jetA++) {
100  _bj_set_NN_crosscheck(jetA, head, jetA);
101  }
102  double * diJ = new double[n];
103  jetA = head;
104  for (int i = 0; i < n; i++) {
105  diJ[i] = _bj_diJ(jetA);
106  jetA++; // have jetA follow i
107  }
108  int history_location = n-1;
109  while (tail != head) {
110  double diJ_min = diJ[0];
111  int diJ_min_jet = 0;
112  for (int i = 1; i < n; i++) {
113  if (diJ[i] < diJ_min) {diJ_min_jet = i; diJ_min = diJ[i];}
114  }
115  history_location++;
116  jetA = & briefjets[diJ_min_jet];
117  jetB = static_cast<BJ *>(jetA->NN);
118  diJ_min *= _invR2;
119  if (jetB != NULL) {
120  if (jetA < jetB) {std::swap(jetA,jetB);}
121  int nn; // new jet index
122  _do_ij_recombination_step(jetA->_jets_index, jetB->_jets_index, diJ_min, nn);
123  _bj_set_jetinfo(jetB, nn);
124  } else {
125  _do_iB_recombination_step(jetA->_jets_index, diJ_min);
126  }
127  tail--; n--;
128  *jetA = *tail;
129  diJ[jetA - head] = diJ[tail-head];
130  for (BJ * jetI = head; jetI != tail; jetI++) {
131  if (jetI->NN == jetA || jetI->NN == jetB) {
132  _bj_set_NN_nocross(jetI, head, tail);
133  diJ[jetI-head] = _bj_diJ(jetI); // update diJ
134  }
135  if (jetB != NULL) {
136  double dist = _bj_dist(jetI,jetB);
137  if (dist < jetI->NN_dist) {
138  if (jetI != jetB) {
139  jetI->NN_dist = dist;
140  jetI->NN = jetB;
141  diJ[jetI-head] = _bj_diJ(jetI); // update diJ...
142  }
143  }
144  if (dist < jetB->NN_dist) {
145  if (jetI != jetB) {
146  jetB->NN_dist = dist;
147  jetB->NN = jetI;}
148  }
149  }
150  if (jetI->NN == tail) {jetI->NN = jetA;}
151  }
152  if (jetB != NULL) {diJ[jetB-head] = _bj_diJ(jetB);}
153  }
154  delete[] diJ;
155  delete[] briefjets;
156 }
158 #endif // __FJCORE_CLUSTERQUENCE_N2_ICC__
159 #ifndef __FJCORE_DYNAMICNEARESTNEIGHBOURS_HH__
160 #define __FJCORE_DYNAMICNEARESTNEIGHBOURS_HH__
161 #include<vector>
162 #include<string>
163 #include<iostream>
164 #include<sstream>
165 #include<cassert>
166 FJCORE_BEGIN_NAMESPACE // defined in fastjet/internal/base.hh
167 class EtaPhi {
168 public:
169  double first, second;
170  EtaPhi() {}
171  EtaPhi(double a, double b) {first = a; second = b;}
172  void sanitize() {
173  if (second < 0) second += twopi;
174  if (second >= twopi) second -= twopi;
175  }
176 };
177 class DnnError : public Error {
178 public:
179  DnnError(const std::string & message_in) : Error(message_in) {}
180 };
182 public:
183  virtual int NearestNeighbourIndex(const int ii) const = 0;
184  virtual double NearestNeighbourDistance(const int ii) const = 0;
185  virtual bool Valid(const int index) const = 0;
186  virtual void RemoveAndAddPoints(const std::vector<int> & indices_to_remove,
187  const std::vector<EtaPhi> & points_to_add,
188  std::vector<int> & indices_added,
189  std::vector<int> & indices_of_updated_neighbours) = 0;
190  inline void RemovePoint (const int index,
191  std::vector<int> & indices_of_updated_neighbours) {
192  std::vector<int> indices_added;
193  std::vector<EtaPhi> points_to_add;
194  std::vector<int> indices_to_remove(1);
195  indices_to_remove[0] = index;
196  RemoveAndAddPoints(indices_to_remove, points_to_add, indices_added,
197  indices_of_updated_neighbours
198  );};
200  const int index1, const int index2,
201  const EtaPhi & newpoint,
202  int & index3,
203  std::vector<int> & indices_of_updated_neighbours) {
204  std::vector<int> indices_added(1);
205  std::vector<EtaPhi> points_to_add(1);
206  std::vector<int> indices_to_remove(2);
207  indices_to_remove[0] = index1;
208  indices_to_remove[1] = index2;
209  points_to_add[0] = newpoint;
210  RemoveAndAddPoints(indices_to_remove, points_to_add, indices_added,
211  indices_of_updated_neighbours
212  );
213  index3 = indices_added[0];
214  };
216 };
218 #endif // __FJCORE_DYNAMICNEARESTNEIGHBOURS_HH__
219 #ifndef __FJCORE_SEARCHTREE_HH__
220 #define __FJCORE_SEARCHTREE_HH__
221 #include<vector>
222 #include<cassert>
223 #include<cstddef>
224 FJCORE_BEGIN_NAMESPACE // defined in fastjet/internal/base.hh
225 template<class T> class SearchTree {
226 public:
227  class Node;
228  class circulator;
229  class const_circulator;
230  SearchTree(const std::vector<T> & init);
231  SearchTree(const std::vector<T> & init, unsigned int max_size);
232  void remove(unsigned node_index);
233  void remove(typename SearchTree::Node * node);
234  void remove(typename SearchTree::circulator & circ);
235  circulator insert(const T & value);
236  const Node & operator[](int i) const {return _nodes[i];};
237  unsigned int size() const {return _nodes.size() - _available_nodes.size();}
238  void verify_structure();
239  void verify_structure_linear() const;
240  void verify_structure_recursive(const Node * , const Node * , const Node * ) const;
241  void print_elements();
242 #ifdef __FJCORE_SEARCHTREE_TRACK_DEPTH
243  inline unsigned int max_depth() const {return _max_depth;};
244 #else
245  inline unsigned int max_depth() const {return 0;};
246 #endif
247  int loc(const Node * node) const ;
248  Node * _find_predecessor(const Node *);
249  Node * _find_successor(const Node *);
250  const Node & operator[](unsigned int i) const {return _nodes[i];};
251  const_circulator somewhere() const;
252  circulator somewhere();
253 private:
254  void _initialize(const std::vector<T> & init);
255  std::vector<Node> _nodes;
256  std::vector<Node *> _available_nodes;
258  unsigned int _n_removes;
259  void _do_initial_connections(unsigned int this_one, unsigned int scale,
260  unsigned int left_edge, unsigned int right_edge,
261  unsigned int depth);
262 #ifdef __FJCORE_SEARCHTREE_TRACK_DEPTH
263  unsigned int _max_depth;
264 #endif
265 };
266 template<class T> class SearchTree<T>::Node{
267 public:
268  Node() {};
269  bool treelinks_null() const {
270  return ((parent==0) && (left==0) && (right==0));};
271  inline void nullify_treelinks() {
272  parent = NULL;
273  left = NULL;
274  right = NULL;
275  };
276  void reset_parents_link_to_me(Node * XX);
283 };
285  if (parent == NULL) {return;}
286  if (parent->right == this) {parent->right = XX;}
287  else {parent->left = XX;}
288 }
289 template<class T> class SearchTree<T>::circulator{
290 public:
291  friend class SearchTree<T>::const_circulator;
292  friend class SearchTree<T>;
293  circulator() : _node(NULL) {}
294  circulator(Node * node) : _node(node) {}
295  const T * operator->() const {return &(_node->value);}
296  T * operator->() {return &(_node->value);}
297  const T & operator*() const {return _node->value;}
298  T & operator*() {return _node->value;}
300  _node = _node->successor;
301  return *this;}
303  circulator tmp = *this;
304  _node = _node->successor;
305  return tmp;}
307  _node = _node->predecessor;
308  return *this;}
310  circulator tmp = *this;
311  _node = _node->predecessor;
312  return tmp;}
313  circulator next() const {
314  return circulator(_node->successor);}
316  return circulator(_node->predecessor);}
317  bool operator!=(const circulator & other) const {return other._node != _node;}
318  bool operator==(const circulator & other) const {return other._node == _node;}
319 private:
321 };
322 template<class T> class SearchTree<T>::const_circulator{
323 public:
324  const_circulator() : _node(NULL) {}
325  const_circulator(const Node * node) : _node(node) {}
326  const_circulator(const circulator & circ) :_node(circ._node) {}
327  const T * operator->() {return &(_node->value);}
328  const T & operator*() const {return _node->value;}
330  _node = _node->successor;
331  return *this;}
333  const_circulator tmp = *this;
334  _node = _node->successor;
335  return tmp;}
337  _node = _node->predecessor;
338  return *this;}
340  const_circulator tmp = *this;
341  _node = _node->predecessor;
342  return tmp;}
344  return const_circulator(_node->successor);}
346  return const_circulator(_node->predecessor);}
347  bool operator!=(const const_circulator & other) const {return other._node != _node;}
348  bool operator==(const const_circulator & other) const {return other._node == _node;}
349 private:
350  const Node * _node;
351 };
352 template<class T> SearchTree<T>::SearchTree(const std::vector<T> & init,
353  unsigned int max_size) :
354  _nodes(max_size) {
355  _available_nodes.reserve(max_size);
356  _available_nodes.resize(max_size - init.size());
357  for (unsigned int i = init.size(); i < max_size; i++) {
358  _available_nodes[i-init.size()] = &(_nodes[i]);
359  }
360  _initialize(init);
361 }
362 template<class T> SearchTree<T>::SearchTree(const std::vector<T> & init) :
363  _nodes(init.size()), _available_nodes(0) {
364  _available_nodes.reserve(init.size());
365  _initialize(init);
366 }
367 template<class T> void SearchTree<T>::_initialize(const std::vector<T> & init) {
368  _n_removes = 0;
369  unsigned n = init.size();
370  assert(n>=1);
371 #ifdef __FJCORE_SEARCHTREE_TRACK_DEPTH
372  _max_depth = 0;
373 #endif
374  for (unsigned int i = 1; i<n; i++) {
375  assert(!(init[i] < init[i-1]));
376  }
377  for(unsigned int i = 0; i < n; i++) {
378  _nodes[i].value = init[i];
379  _nodes[i].predecessor = (& (_nodes[i])) - 1;
380  _nodes[i].successor = (& (_nodes[i])) + 1;
381  _nodes[i].nullify_treelinks();
382  }
383  _nodes[0].predecessor = (& (_nodes[n-1]));
384  _nodes[n-1].successor = (& (_nodes[0]));
385  unsigned int scale = (n+1)/2;
386  unsigned int top = std::min(n-1,scale);
387  _nodes[top].parent = NULL;
388  _top_node = &(_nodes[top]);
389  _do_initial_connections(top, scale, 0, n, 0);
390 }
391 template<class T> inline int SearchTree<T>::loc(const Node * node) const {return node == NULL?
392  -999 : node - &(_nodes[0]);}
394  unsigned int this_one,
395  unsigned int scale,
396  unsigned int left_edge,
397  unsigned int right_edge,
398  unsigned int depth
399  ) {
400 #ifdef __FJCORE_SEARCHTREE_TRACK_DEPTH
401  _max_depth = max(depth, _max_depth);
402 #endif
403  unsigned int ref_new_scale = (scale+1)/2;
404  unsigned new_scale = ref_new_scale;
405  bool did_child = false;
406  while(true) {
407  int left = this_one - new_scale; // be careful here to use signed int...
408  if (left >= static_cast<int>(left_edge)
409  && _nodes[left].treelinks_null() ) {
410  _nodes[left].parent = &(_nodes[this_one]);
411  _nodes[this_one].left = &(_nodes[left]);
412  _do_initial_connections(left, new_scale, left_edge, this_one, depth+1);
413  did_child = true;
414  break;
415  }
416  unsigned int old_new_scale = new_scale;
417  new_scale = (old_new_scale + 1)/2;
418  if (new_scale == old_new_scale) break;
419  }
420  if (!did_child) {_nodes[this_one].left = NULL;}
421  new_scale = ref_new_scale;
422  did_child = false;
423  while(true) {
424  unsigned int right = this_one + new_scale;
425  if (right < right_edge && _nodes[right].treelinks_null()) {
426  _nodes[right].parent = &(_nodes[this_one]);
427  _nodes[this_one].right = &(_nodes[right]);
428  _do_initial_connections(right, new_scale, this_one+1,right_edge,depth+1);
429  did_child = true;
430  break;
431  }
432  unsigned int old_new_scale = new_scale;
433  new_scale = (old_new_scale + 1)/2;
434  if (new_scale == old_new_scale) break;
435  }
436  if (!did_child) {_nodes[this_one].right = NULL;}
437 }
438 template<class T> void SearchTree<T>::remove(unsigned int node_index) {
439  remove(&(_nodes[node_index]));
440 }
441 template<class T> void SearchTree<T>::remove(circulator & circ) {
442  remove(circ._node);
443 }
444 template<class T> void SearchTree<T>::remove(typename SearchTree<T>::Node * node) {
445  assert(size() > 1); // switch this to throw...?
446  assert(!node->treelinks_null());
447  node->predecessor->successor = node->successor;
448  node->successor->predecessor = node->predecessor;
449  if (node->left == NULL && node->right == NULL) {
450  node->reset_parents_link_to_me(NULL);
451  } else if (node->left != NULL && node->right == NULL){
452  node->reset_parents_link_to_me(node->left);
453  node->left->parent = node->parent;
454  if (_top_node == node) {_top_node = node->left;}
455  } else if (node->left == NULL && node->right != NULL){
456  node->reset_parents_link_to_me(node->right);
457  node->right->parent = node->parent;
458  if (_top_node == node) {_top_node = node->right;}
459  } else {
460  Node * replacement;
461  bool use_predecessor = (_n_removes % 2 == 1);
462  if (use_predecessor) {
463  replacement = node->predecessor;
464  assert(replacement->right == NULL); // guaranteed if it's our predecessor
465  if (replacement != node->left) {
466  if (replacement->left != NULL) {
467  replacement->left->parent = replacement->parent;}
468  replacement->reset_parents_link_to_me(replacement->left);
469  replacement->left = node->left;
470  }
471  replacement->parent = node->parent;
472  replacement->right = node->right;
473  } else {
474  replacement = node->successor;
475  assert(replacement->left == NULL); // guaranteed if it's our successor
476  if (replacement != node->right) {
477  if (replacement->right != NULL) {
478  replacement->right->parent = replacement->parent;}
479  replacement->reset_parents_link_to_me(replacement->right);
480  replacement->right = node->right;
481  }
482  replacement->parent = node->parent;
483  replacement->left = node->left;
484  }
485  node->reset_parents_link_to_me(replacement);
486  if (node->left != replacement) {node->left->parent = replacement;}
487  if (node->right != replacement) {node->right->parent = replacement;}
488  if (_top_node == node) {_top_node = replacement;}
489  }
490  node->nullify_treelinks();
491  node->predecessor = NULL;
492  node->successor = NULL;
493  _n_removes++;
494  _available_nodes.push_back(node);
495 }
496 template<class T> typename SearchTree<T>::circulator SearchTree<T>::insert(const T & value) {
497  assert(_available_nodes.size() > 0);
498  Node * node = _available_nodes.back();
499  _available_nodes.pop_back();
500  node->value = value;
501  Node * location = _top_node;
502  Node * old_location = NULL;
503  bool on_left = true; // (init not needed -- but soothes g++4)
504 #ifdef __FJCORE_SEARCHTREE_TRACK_DEPTH
505  unsigned int depth = 0;
506 #endif
507  while(location != NULL) {
508 #ifdef __FJCORE_SEARCHTREE_TRACK_DEPTH
509  depth++;
510 #endif
511  old_location = location;
512  on_left = value < location->value;
513  if (on_left) {location = location->left;}
514  else {location = location->right;}
515  }
516 #ifdef __FJCORE_SEARCHTREE_TRACK_DEPTH
517  _max_depth = max(depth, _max_depth);
518 #endif
519  node->parent = old_location;
520  if (on_left) {node->parent->left = node;}
521  else {node->parent->right = node;}
522  node->left = NULL;
523  node->right = NULL;
524  node->predecessor = _find_predecessor(node);
525  if (node->predecessor != NULL) {
526  node->successor = node->predecessor->successor;
527  node->predecessor->successor = node;
528  node->successor->predecessor = node;
529  } else {
530  node->successor = _find_successor(node);
531  assert(node->successor != NULL); // can only happen if we're sole element
532  node->predecessor = node->successor->predecessor;
533  node->successor->predecessor = node;
534  node->predecessor->successor = node;
535  }
536  return circulator(node);
537 }
538 template<class T> void SearchTree<T>::verify_structure() {
539  verify_structure_linear();
540  const Node * left_limit = _top_node;
541  while (left_limit->left != NULL) {left_limit = left_limit->left;}
542  const Node * right_limit = _top_node;
543  while (right_limit->right != NULL) {right_limit = right_limit->right;}
544  verify_structure_recursive(_top_node, left_limit, right_limit);
545 }
547  const typename SearchTree<T>::Node * element,
548  const typename SearchTree<T>::Node * left_limit,
549  const typename SearchTree<T>::Node * right_limit) const {
550  assert(!(element->value < left_limit->value));
551  assert(!(right_limit->value < element->value));
552  const Node * left = element->left;
553  if (left != NULL) {
554  assert(!(element->value < left->value));
555  if (left != left_limit) {
556  verify_structure_recursive(left, left_limit, element);}
557  }
558  const Node * right = element->right;
559  if (right != NULL) {
560  assert(!(right->value < element->value));
561  if (right != right_limit) {
562  verify_structure_recursive(right, element, right_limit);}
563  }
564 }
565 template<class T> void SearchTree<T>::verify_structure_linear() const {
566  unsigned n_top = 0;
567  unsigned n_null = 0;
568  for(unsigned i = 0; i < _nodes.size(); i++) {
569  const typename SearchTree<T>::Node * node = &(_nodes[i]);
570  if (node->treelinks_null()) {n_null++; continue;}
571  if (node->parent == NULL) {
572  n_top++;
573  } else {
574  assert((node->parent->left == node) ^ (node->parent->right == node));
575  }
576  if (node->left != NULL) {
577  assert(!(node->value < node->left->value ));}
578  if (node->right != NULL) {
579  assert(!(node->right->value < node->value ));}
580  }
581  assert(n_top == 1 || (n_top == 0 && size() <= 1) );
582  assert(n_null == _available_nodes.size() ||
583  (n_null == _available_nodes.size() + 1 && size() == 1));
584 }
585 template<class T> typename SearchTree<T>::Node * SearchTree<T>::_find_predecessor(const typename SearchTree<T>::Node * node) {
586  typename SearchTree<T>::Node * newnode;
587  if (node->left != NULL) {
588  newnode = node->left;
589  while(newnode->right != NULL) {newnode = newnode->right;}
590  return newnode;
591  } else {
592  const typename SearchTree<T>::Node * lastnode = node;
593  newnode = node->parent;
594  while(newnode != NULL) {
595  if (newnode->right == lastnode) {return newnode;}
596  lastnode = newnode;
597  newnode = newnode->parent;
598  }
599  return newnode;
600  }
601 }
602 template<class T> typename SearchTree<T>::Node * SearchTree<T>::_find_successor(const typename SearchTree<T>::Node * node) {
603  typename SearchTree<T>::Node * newnode;
604  if (node->right != NULL) {
605  newnode = node->right;
606  while(newnode->left != NULL) {newnode = newnode->left;}
607  return newnode;
608  } else {
609  const typename SearchTree<T>::Node * lastnode = node;
610  newnode = node->parent;
611  while(newnode != NULL) {
612  if (newnode->left == lastnode) {return newnode;}
613  lastnode = newnode;
614  newnode = newnode->parent;
615  }
616  return newnode;
617  }
618 }
619 template<class T> void SearchTree<T>::print_elements() {
620  typename SearchTree<T>::Node * base_node = &(_nodes[0]);
621  typename SearchTree<T>::Node * node = base_node;
622  int n = _nodes.size();
623  for(; node - base_node < n ; node++) {
624  printf("%4d parent:%4d left:%4d right:%4d pred:%4d succ:%4d value:%10.6f\n",loc(node), loc(node->parent), loc(node->left), loc(node->right), loc(node->predecessor),loc(node->successor),node->value);
625  }
626 }
628  return circulator(_top_node);
629 }
630 template<class T> typename SearchTree<T>::const_circulator SearchTree<T>::somewhere() const {
631  return const_circulator(_top_node);
632 }
634 #endif // __FJCORE_SEARCHTREE_HH__
635 #ifndef __FJCORE_MINHEAP__HH__
636 #define __FJCORE_MINHEAP__HH__
637 #include<vector>
638 #include<cassert>
639 #include<memory>
640 #include<limits>
641 FJCORE_BEGIN_NAMESPACE // defined in fastjet/internal/base.hh
642 class MinHeap {
643 public:
644  MinHeap (const std::vector<double> & values, unsigned int max_size) :
645  _heap(max_size) {initialise(values);}
646  MinHeap (unsigned int max_size) : _heap(max_size) {}
647  MinHeap (const std::vector<double> & values) :
648  _heap(values.size()) {initialise(values);}
649  void initialise(const std::vector<double> & values);
650  inline unsigned int minloc() const {
651  return (_heap[0].minloc) - &(_heap[0]);}
652  inline double minval() const {return _heap[0].minloc->value;}
653  inline double operator[](int i) const {return _heap[i].value;}
654  void remove(unsigned int loc) {
655  update(loc,std::numeric_limits<double>::max());};
656  void update(unsigned int, double);
657 private:
658  struct ValueLoc{
659  double value;
661  };
662  std::vector<ValueLoc> _heap;
663 };
665 #endif // __FJCORE_MINHEAP__HH__
666 #ifndef __FJCORE_CLOSESTPAIR2DBASE__HH__
667 #define __FJCORE_CLOSESTPAIR2DBASE__HH__
668 #include<vector>
669 FJCORE_BEGIN_NAMESPACE // defined in fastjet/internal/base.hh
670 class Coord2D {
671 public:
672  double x, y;
673  Coord2D() : x(0.0), y(0.0) {};
674  Coord2D(double a, double b): x(a), y(b) {};
675  Coord2D operator-(const Coord2D & other) const {
676  return Coord2D(x - other.x, y - other.y);};
677  Coord2D operator+(const Coord2D & other) const {
678  return Coord2D(x + other.x, y + other.y);};
679  Coord2D operator*(double factor) const {return Coord2D(factor*x,factor*y);};
680  friend Coord2D operator*(double factor, const Coord2D & coord) {
681  return Coord2D(factor*coord.x,factor*coord.y);
682  }
683  Coord2D operator/(double divisor) const {
684  return Coord2D(x / divisor, y / divisor);};
685  friend double distance2(const Coord2D & a, const Coord2D & b) {
686  double dx = a.x - b.x, dy = a.y-b.y;
687  return dx*dx+dy*dy;
688  };
689  double distance2(const Coord2D & b) const {
690  double dx = x - b.x, dy = y-b.y;
691  return dx*dx+dy*dy;
692  };
693 };
695 public:
696  virtual void closest_pair(unsigned int & ID1, unsigned int & ID2,
697  double & distance2) const = 0;
698  virtual void remove(unsigned int ID) = 0;
699  virtual unsigned int insert(const Coord2D & position) = 0;
700  virtual unsigned int replace(unsigned int ID1, unsigned int ID2,
701  const Coord2D & position) {
702  remove(ID1);
703  remove(ID2);
704  unsigned new_ID = insert(position);
705  return(new_ID);
706  };
707  virtual void replace_many(const std::vector<unsigned int> & IDs_to_remove,
708  const std::vector<Coord2D> & new_positions,
709  std::vector<unsigned int> & new_IDs) {
710  for(unsigned i = 0; i < IDs_to_remove.size(); i++) {
711  remove(IDs_to_remove[i]);}
712  new_IDs.resize(0);
713  for(unsigned i = 0; i < new_positions.size(); i++) {
714  new_IDs.push_back(insert(new_positions[i]));}
715  }
716  virtual unsigned int size() = 0;
717  virtual ~ClosestPair2DBase() {};
718 };
720 #endif // __FJCORE_CLOSESTPAIR2DBASE__HH__
721 #ifndef __FJCORE_CLOSESTPAIR2D__HH__
722 #define __FJCORE_CLOSESTPAIR2D__HH__
723 #include<vector>
724 #include<stack>
725 #include<iostream>
726 FJCORE_BEGIN_NAMESPACE // defined in fastjet/internal/base.hh
728 public:
729  ClosestPair2D(const std::vector<Coord2D> & positions,
730  const Coord2D & left_corner, const Coord2D & right_corner) {
731  _initialize(positions, left_corner, right_corner, positions.size());
732  };
733  ClosestPair2D(const std::vector<Coord2D> & positions,
734  const Coord2D & left_corner, const Coord2D & right_corner,
735  const unsigned int max_size) {
736  _initialize(positions, left_corner, right_corner, max_size);
737  };
738  void closest_pair(unsigned int & ID1, unsigned int & ID2,
739  double & distance2) const;
740  void remove(unsigned int ID);
741  unsigned int insert(const Coord2D &);
742  virtual unsigned int replace(unsigned int ID1, unsigned int ID2,
743  const Coord2D & position);
744  virtual void replace_many(const std::vector<unsigned int> & IDs_to_remove,
745  const std::vector<Coord2D> & new_positions,
746  std::vector<unsigned int> & new_IDs);
747  inline void print_tree_depths(std::ostream & outdev) const {
748  outdev << _trees[0]->max_depth() << " "
749  << _trees[1]->max_depth() << " "
750  << _trees[2]->max_depth() << "\n";
751  };
752  unsigned int size();
753 private:
754  void _initialize(const std::vector<Coord2D> & positions,
755  const Coord2D & left_corner, const Coord2D & right_corner,
756  const unsigned int max_size);
757  static const unsigned int _nshift = 3;
758  class Point; // will be defined below
759  template<class T> class triplet {
760  public:
761  inline const T & operator[](unsigned int i) const {return _contents[i];};
762  inline T & operator[](unsigned int i) {return _contents[i];};
763  private:
765  };
766  class Shuffle {
767  public:
768  unsigned int x, y;
770  bool operator<(const Shuffle &) const;
771  void operator+=(unsigned int shift) {x += shift; y+= shift;};
772  };
778  std::vector<Point> _points;
779  std::stack<Point *> _available_points;
780  std::vector<Point *> _points_under_review;
781  static const unsigned int _remove_heap_entry = 1;
782  static const unsigned int _review_heap_entry = 2;
783  static const unsigned int _review_neighbour = 4;
784  void _add_label(Point * point, unsigned int review_flag);
785  void _set_label(Point * point, unsigned int review_flag);
787  void _remove_from_search_tree(Point * point_to_remove);
788  void _insert_into_search_tree(Point * new_point);
789  void _point2shuffle(Point & , Shuffle & , unsigned int shift);
791  double _range;
792  int _ID(const Point *) const;
793  triplet<unsigned int> _shifts; // absolute shifts
794  triplet<unsigned int> _rel_shifts; // shifts relative to previous shift
795  unsigned int _cp_search_range;
796 };
798 public:
803  unsigned int review_flag;
804  double distance2(const Point & other) const {
805  return coord.distance2(other.coord);
806  };
807 };
808 inline bool floor_ln2_less(unsigned x, unsigned y) {
809  if (x>y) return false;
810  return (x < (x^y)); // beware of operator precedence...
811 }
812 inline int ClosestPair2D::_ID(const Point * point) const {
813  return point - &(_points[0]);
814 }
815 inline unsigned int ClosestPair2D::size() {
816  return _points.size() - _available_points.size();
817 }
819 #endif // __FJCORE_CLOSESTPAIR2D__HH__
820 #ifndef __FJCORE_LAZYTILING9ALT_HH__
821 #define __FJCORE_LAZYTILING9ALT_HH__
822 FJCORE_BEGIN_NAMESPACE // defined in fastjet/internal/base.hh
823 const double tile_edge_security_margin=1.0e-7;
824 class TiledJet {
825 public:
826  double eta, phi, kt2, NN_dist;
832  inline bool minheap_update_needed() const {return _minheap_update_needed;}
833 };
834 const int n_tile_neighbours = 9;
835 class Tile {
836 public:
837  typedef double (Tile::*DistToTileFn)(const TiledJet*) const;
838  typedef std::pair<Tile *, DistToTileFn> TileFnPair;
844  bool tagged;
846  double max_NN_dist;
848  double distance_to_centre(const TiledJet *) const {return 0;}
849  double distance_to_left(const TiledJet * jet) const {
850  double deta = jet->eta - eta_min;
851  return deta*deta;
852  }
853  double distance_to_right(const TiledJet * jet) const {
854  double deta = jet->eta - eta_max;
855  return deta*deta;
856  }
857  double distance_to_bottom(const TiledJet * jet) const {
858  double dphi = jet->phi - phi_min;
859  return dphi*dphi;
860  }
861  double distance_to_top(const TiledJet * jet) const {
862  double dphi = jet->phi - phi_max;
863  return dphi*dphi;
864  }
865  double distance_to_left_top(const TiledJet * jet) const {
866  double deta = jet->eta - eta_min;
867  double dphi = jet->phi - phi_max;
868  return deta*deta + dphi*dphi;
869  }
870  double distance_to_left_bottom(const TiledJet * jet) const {
871  double deta = jet->eta - eta_min;
872  double dphi = jet->phi - phi_min;
873  return deta*deta + dphi*dphi;
874  }
875  double distance_to_right_top(const TiledJet * jet) const {
876  double deta = jet->eta - eta_max;
877  double dphi = jet->phi - phi_max;
878  return deta*deta + dphi*dphi;
879  }
880  double distance_to_right_bottom(const TiledJet * jet) const {
881  double deta = jet->eta - eta_max;
882  double dphi = jet->phi - phi_min;
883  return deta*deta + dphi*dphi;
884  }
885 };
887 public:
888  LazyTiling9Alt(ClusterSequence & cs);
889  void run();
890 protected:
891  ClusterSequence & _cs;
892  const std::vector<PseudoJet> & _jets;
893  std::vector<Tile> _tiles;
894  double _Rparam, _R2, _invR2;
899  std::vector<TiledJet *> _jets_for_minheap;
900  void _initialise_tiles();
901  inline int _tile_index (int ieta, int iphi) const {
902  return (ieta-_tiles_ieta_min)*_n_tiles_phi
903  + (iphi+_n_tiles_phi) % _n_tiles_phi;
904  }
905  void _bj_remove_from_tiles(TiledJet * const jet);
906  int _tile_index(const double eta, const double phi) const;
907  void _tj_set_jetinfo(TiledJet * const jet, const int _jets_index);
908  void _print_tiles(TiledJet * briefjets ) const;
909  void _add_neighbours_to_tile_union(const int tile_index,
910  std::vector<int> & tile_union, int & n_near_tiles) const;
911  void _add_untagged_neighbours_to_tile_union(const int tile_index,
912  std::vector<int> & tile_union, int & n_near_tiles);
914  std::vector<int> & tile_union, int & n_near_tiles);
915  void _update_jetX_jetI_NN(TiledJet * jetX, TiledJet * jetI, std::vector<TiledJet *> & jets_for_minheap);
916  void _set_NN(TiledJet * jetI, std::vector<TiledJet *> & jets_for_minheap);
917  template <class J> double _bj_diJ(const J * const jet) const {
918  double kt2 = jet->kt2;
919  if (jet->NN != NULL) {if (jet->NN->kt2 < kt2) {kt2 = jet->NN->kt2;}}
920  return jet->NN_dist * kt2;
921  }
922  template <class J> inline void _bj_set_jetinfo(
923  J * const jetA, const int _jets_index) const {
924  jetA->eta = _jets[_jets_index].rap();
925  jetA->phi = _jets[_jets_index].phi_02pi();
926  jetA->kt2 = _cs.jet_scale_for_algorithm(_jets[_jets_index]);
927  jetA->_jets_index = _jets_index;
928  jetA->NN_dist = _R2;
929  jetA->NN = NULL;
930  }
931  template <class J> inline double _bj_dist(
932  const J * const jetA, const J * const jetB) const {
933  double dphi = std::abs(jetA->phi - jetB->phi);
934  double deta = (jetA->eta - jetB->eta);
935  if (dphi > pi) {dphi = twopi - dphi;}
936  return dphi*dphi + deta*deta;
937  }
938  template <class J> inline double _bj_dist_not_periodic(
939  const J * const jetA, const J * const jetB) const {
940  double dphi = jetA->phi - jetB->phi;
941  double deta = (jetA->eta - jetB->eta);
942  return dphi*dphi + deta*deta;
943  }
944 };
946 #endif // __FJCORE_LAZYTILING9ALT_HH__
947 #ifndef __FJCORE_LAZYTILING9_HH__
948 #define __FJCORE_LAZYTILING9_HH__
949 FJCORE_BEGIN_NAMESPACE // defined in fastjet/internal/base.hh
950 template<int NN>
951 class Tile2Base {
952 public:
958  bool tagged;
960  double max_NN_dist;
962  int jet_count() const {
963  int count = 0;
964  const TiledJet * jet = head;
965  while (jet != 0) {
966  count++;
967  jet = jet->next;
968  }
969  return count;
970  }
971 };
973 class LazyTiling9 {
974 public:
975  LazyTiling9(ClusterSequence & cs);
976  void run();
977 protected:
978  ClusterSequence & _cs;
979  const std::vector<PseudoJet> & _jets;
980  std::vector<Tile2> _tiles;
981 #ifdef INSTRUMENT2
982  int _ncall; // GPS tmp
983  int _ncall_dtt; // GPS tmp
984 #endif // INSTRUMENT2
985  double _Rparam, _R2, _invR2;
990  std::vector<TiledJet *> _jets_for_minheap;
991  void _initialise_tiles();
992  inline int _tile_index (int ieta, int iphi) const {
993  return (ieta-_tiles_ieta_min)*_n_tiles_phi
994  + (iphi+_n_tiles_phi) % _n_tiles_phi;
995  }
996  void _bj_remove_from_tiles(TiledJet * const jet);
997  int _tile_index(const double eta, const double phi) const;
998  void _tj_set_jetinfo(TiledJet * const jet, const int _jets_index);
999  void _print_tiles(TiledJet * briefjets ) const;
1000  void _add_neighbours_to_tile_union(const int tile_index,
1001  std::vector<int> & tile_union, int & n_near_tiles) const;
1002  void _add_untagged_neighbours_to_tile_union(const int tile_index,
1003  std::vector<int> & tile_union, int & n_near_tiles);
1005  std::vector<int> & tile_union, int & n_near_tiles);
1006  double _distance_to_tile(const TiledJet * bj, const Tile2 *)
1007 #ifdef INSTRUMENT2
1008  ;
1009 #else
1010  const;
1011 #endif
1012  void _update_jetX_jetI_NN(TiledJet * jetX, TiledJet * jetI, std::vector<TiledJet *> & jets_for_minheap);
1013  void _set_NN(TiledJet * jetI, std::vector<TiledJet *> & jets_for_minheap);
1014  template <class J> double _bj_diJ(const J * const jet) const {
1015  double kt2 = jet->kt2;
1016  if (jet->NN != NULL) {if (jet->NN->kt2 < kt2) {kt2 = jet->NN->kt2;}}
1017  return jet->NN_dist * kt2;
1018  }
1019  template <class J> inline void _bj_set_jetinfo(
1020  J * const jetA, const int _jets_index) const {
1021  jetA->eta = _jets[_jets_index].rap();
1022  jetA->phi = _jets[_jets_index].phi_02pi();
1023  jetA->kt2 = _cs.jet_scale_for_algorithm(_jets[_jets_index]);
1024  jetA->_jets_index = _jets_index;
1025  jetA->NN_dist = _R2;
1026  jetA->NN = NULL;
1027  }
1028  template <class J> inline double _bj_dist(
1029  const J * const jetA, const J * const jetB)
1030 #ifdef INSTRUMENT2
1031  {
1032  _ncall++; // GPS tmp
1033 #else
1034  const {
1035 #endif
1036  double dphi = std::abs(jetA->phi - jetB->phi);
1037  double deta = (jetA->eta - jetB->eta);
1038  if (dphi > pi) {dphi = twopi - dphi;}
1039  return dphi*dphi + deta*deta;
1040  }
1041  template <class J> inline double _bj_dist_not_periodic(
1042  const J * const jetA, const J * const jetB)
1043 #ifdef INSTRUMENT2
1044  {
1045  _ncall++; // GPS tmp
1046 #else
1047  const {
1048 #endif
1049  double dphi = jetA->phi - jetB->phi;
1050  double deta = (jetA->eta - jetB->eta);
1051  return dphi*dphi + deta*deta;
1052  }
1053 };
1055 #endif // __FJCORE_LAZYTILING9_HH__
1056 #ifndef __FJCORE_LAZYTILING25_HH__
1057 #define __FJCORE_LAZYTILING25_HH__
1058 FJCORE_BEGIN_NAMESPACE // defined in fastjet/internal/base.hh
1061 public:
1062  LazyTiling25(ClusterSequence & cs);
1063  void run();
1064 protected:
1065  ClusterSequence & _cs;
1066  const std::vector<PseudoJet> & _jets;
1067  std::vector<Tile25> _tiles;
1068 #ifdef INSTRUMENT2
1069  int _ncall; // GPS tmp
1070  int _ncall_dtt; // GPS tmp
1071 #endif // INSTRUMENT2
1072  double _Rparam, _R2, _invR2;
1077  std::vector<TiledJet *> _jets_for_minheap;
1078  void _initialise_tiles();
1079  inline int _tile_index (int ieta, int iphi) const {
1080  return (ieta-_tiles_ieta_min)*_n_tiles_phi
1081  + (iphi+_n_tiles_phi) % _n_tiles_phi;
1082  }
1083  void _bj_remove_from_tiles(TiledJet * const jet);
1084  int _tile_index(const double eta, const double phi) const;
1085  void _tj_set_jetinfo(TiledJet * const jet, const int _jets_index);
1086  void _print_tiles(TiledJet * briefjets ) const;
1087  void _add_neighbours_to_tile_union(const int tile_index,
1088  std::vector<int> & tile_union, int & n_near_tiles) const;
1089  void _add_untagged_neighbours_to_tile_union(const int tile_index,
1090  std::vector<int> & tile_union, int & n_near_tiles);
1092  std::vector<int> & tile_union, int & n_near_tiles);
1093  double _distance_to_tile(const TiledJet * bj, const Tile25 *)
1094 #ifdef INSTRUMENT2
1095  ;
1096 #else
1097  const;
1098 #endif
1099  void _update_jetX_jetI_NN(TiledJet * jetX, TiledJet * jetI, std::vector<TiledJet *> & jets_for_minheap);
1100  void _set_NN(TiledJet * jetI, std::vector<TiledJet *> & jets_for_minheap);
1101  template <class J> double _bj_diJ(const J * const jet) const {
1102  double kt2 = jet->kt2;
1103  if (jet->NN != NULL) {if (jet->NN->kt2 < kt2) {kt2 = jet->NN->kt2;}}
1104  return jet->NN_dist * kt2;
1105  }
1106  template <class J> inline void _bj_set_jetinfo(
1107  J * const jetA, const int _jets_index) const {
1108  jetA->eta = _jets[_jets_index].rap();
1109  jetA->phi = _jets[_jets_index].phi_02pi();
1110  jetA->kt2 = _cs.jet_scale_for_algorithm(_jets[_jets_index]);
1111  jetA->_jets_index = _jets_index;
1112  jetA->NN_dist = _R2;
1113  jetA->NN = NULL;
1114  }
1115  template <class J> inline double _bj_dist(
1116  const J * const jetA, const J * const jetB)
1117 #ifdef INSTRUMENT2
1118  {
1119  _ncall++; // GPS tmp
1120 #else
1121  const {
1122 #endif
1123  double dphi = std::abs(jetA->phi - jetB->phi);
1124  double deta = (jetA->eta - jetB->eta);
1125  if (dphi > pi) {dphi = twopi - dphi;}
1126  return dphi*dphi + deta*deta;
1127  }
1128  template <class J> inline double _bj_dist_not_periodic(
1129  const J * const jetA, const J * const jetB)
1130 #ifdef INSTRUMENT2
1131  {
1132  _ncall++; // GPS tmp
1133 #else
1134  const {
1135 #endif
1136  double dphi = jetA->phi - jetB->phi;
1137  double deta = (jetA->eta - jetB->eta);
1138  return dphi*dphi + deta*deta;
1139  }
1140 };
1142 #endif // __FJCORE_LAZYTILING25_HH__
1143 #ifndef __FJCORE_TILINGEXTENT_HH__
1144 #define __FJCORE_TILINGEXTENT_HH__
1145 FJCORE_BEGIN_NAMESPACE // defined in fastjet/internal/base.hh
1147 public:
1148  TilingExtent(ClusterSequence & cs);
1149  TilingExtent(const std::vector<PseudoJet> &particles);
1150  double minrap() const {return _minrap;}
1151  double maxrap() const {return _maxrap;}
1152  double sum_of_binned_squared_multiplicity() const {return _cumul2;}
1153 private:
1154  double _minrap, _maxrap, _cumul2;
1155  void _determine_rapidity_extent(const std::vector<PseudoJet> & particles);
1156 };
1157 FJCORE_END_NAMESPACE // defined in fastjet/internal/base.hh
1158 #endif // __FJCORE_TILINGEXTENT_HH__
1159 #include<limits>
1160 #include<iostream>
1161 #include<iomanip>
1162 #include<algorithm>
1163 FJCORE_BEGIN_NAMESPACE // defined in fastjet/internal/base.hh
1164 const unsigned int twopow31 = 2147483648U;
1165 using namespace std;
1167  unsigned int shift) {
1168  Coord2D renorm_point = (point.coord - _left_corner)/_range;
1169  assert(renorm_point.x >=0);
1170  assert(renorm_point.x <=1);
1171  assert(renorm_point.y >=0);
1172  assert(renorm_point.y <=1);
1173  shuffle.x = static_cast<unsigned int>(twopow31 * renorm_point.x) + shift;
1174  shuffle.y = static_cast<unsigned int>(twopow31 * renorm_point.y) + shift;
1175  shuffle.point = &point;
1176 }
1178  if (floor_ln2_less(x ^ q.x, y ^ q.y)) {
1179  return (y < q.y);
1180  } else {
1181  return (x < q.x);
1182  }
1183 }
1184 void ClosestPair2D::_initialize(const std::vector<Coord2D> & positions,
1185  const Coord2D & left_corner,
1186  const Coord2D & right_corner,
1187  unsigned int max_size) {
1188  unsigned int n_positions = positions.size();
1189  assert(max_size >= n_positions);
1190  _points.resize(max_size);
1191  for (unsigned int i = n_positions; i < max_size; i++) {
1192  _available_points.push(&(_points[i]));
1193  }
1194  _left_corner = left_corner;
1195  _range = max((right_corner.x - left_corner.x),
1196  (right_corner.y - left_corner.y));
1197  vector<Shuffle> shuffles(n_positions);
1198  for (unsigned int i = 0; i < n_positions; i++) {
1199  _points[i].coord = positions[i];
1200  _points[i].neighbour_dist2 = numeric_limits<double>::max();
1201  _points[i].review_flag = 0;
1202  _point2shuffle(_points[i], shuffles[i], 0);
1203  }
1204  for (unsigned ishift = 0; ishift < _nshift; ishift++) {
1205  _shifts[ishift] = static_cast<unsigned int>(((twopow31*1.0)*ishift)/_nshift);
1206  if (ishift == 0) {_rel_shifts[ishift] = 0;}
1207  else {_rel_shifts[ishift] = _shifts[ishift] - _shifts[ishift-1];}
1208  }
1209  _cp_search_range = 30;
1210  _points_under_review.reserve(_nshift * _cp_search_range);
1211  for (unsigned int ishift = 0; ishift < _nshift; ishift++) {
1212  if (ishift > 0) {
1213  unsigned rel_shift = _rel_shifts[ishift];
1214  for (unsigned int i = 0; i < shuffles.size(); i++) {
1215  shuffles[i] += rel_shift; }
1216  }
1217  sort(shuffles.begin(), shuffles.end());
1218  _trees[ishift] = SharedPtr<Tree>(new Tree(shuffles, max_size));
1219  circulator circ = _trees[ishift]->somewhere(), start=circ;
1220  unsigned int CP_range = min(_cp_search_range, n_positions-1);
1221  do {
1222  Point * this_point = circ->point;
1223  this_point->circ[ishift] = circ;
1224  circulator other = circ;
1225  for (unsigned i=0; i < CP_range; i++) {
1226  ++other;
1227  double dist2 = this_point->distance2(*other->point);
1228  if (dist2 < this_point->neighbour_dist2) {
1229  this_point->neighbour_dist2 = dist2;
1230  this_point->neighbour = other->point;
1231  }
1232  }
1233  } while (++circ != start);
1234  }
1235  vector<double> mindists2(n_positions);
1236  for (unsigned int i = 0; i < n_positions; i++) {
1237  mindists2[i] = _points[i].neighbour_dist2;}
1238  _heap = SharedPtr<MinHeap>(new MinHeap(mindists2, max_size));
1239 }
1240 void ClosestPair2D::closest_pair(unsigned int & ID1, unsigned int & ID2,
1241  double & distance2) const {
1242  ID1 = _heap->minloc();
1243  ID2 = _ID(_points[ID1].neighbour);
1244  distance2 = _points[ID1].neighbour_dist2;
1245  if (ID1 > ID2) std::swap(ID1,ID2);
1246 }
1247 inline void ClosestPair2D::_add_label(Point * point, unsigned int review_flag) {
1248  if (point->review_flag == 0) _points_under_review.push_back(point);
1249  point->review_flag |= review_flag;
1250 }
1251 inline void ClosestPair2D::_set_label(Point * point, unsigned int review_flag) {
1252  if (point->review_flag == 0) _points_under_review.push_back(point);
1253  point->review_flag = review_flag;
1254 }
1255 void ClosestPair2D::remove(unsigned int ID) {
1256  Point * point_to_remove = & (_points[ID]);
1257  _remove_from_search_tree(point_to_remove);
1258  _deal_with_points_to_review();
1259 }
1261  _available_points.push(point_to_remove);
1262  _set_label(point_to_remove, _remove_heap_entry);
1263  unsigned int CP_range = min(_cp_search_range, size()-1);
1264  for (unsigned int ishift = 0; ishift < _nshift; ishift++) {
1265  circulator removed_circ = point_to_remove->circ[ishift];
1266  circulator right_end = removed_circ.next();
1267  _trees[ishift]->remove(removed_circ);
1268  circulator left_end = right_end, orig_right_end = right_end;
1269  for (unsigned int i = 0; i < CP_range; i++) {left_end--;}
1270  if (size()-1 < _cp_search_range) {
1271  left_end--; right_end--;
1272  }
1273  do {
1274  Point * left_point = left_end->point;
1275  if (left_point->neighbour == point_to_remove) {
1276  // we'll deal with it later...
1277  _add_label(left_point, _review_neighbour);
1278  } else {
1279  // check to see if right point has become its closest neighbour
1280  double dist2 = left_point->distance2(*right_end->point);
1281  if (dist2 < left_point->neighbour_dist2) {
1282  left_point->neighbour = right_end->point;
1283  left_point->neighbour_dist2 = dist2;
1284  // NB: (LESSER) REVIEW NEEDED HERE TOO...
1285  _add_label(left_point, _review_heap_entry);
1286  }
1287  }
1288  ++right_end;
1289  } while (++left_end != orig_right_end);
1290  } // ishift...
1291 }
1293  unsigned int CP_range = min(_cp_search_range, size()-1);
1294  while(_points_under_review.size() > 0) {
1295  Point * this_point = _points_under_review.back();
1296  _points_under_review.pop_back();
1297  if (this_point->review_flag & _remove_heap_entry) {
1298  assert(!(this_point->review_flag ^ _remove_heap_entry));
1299  _heap->remove(_ID(this_point));
1300  }
1301  else {
1302  if (this_point->review_flag & _review_neighbour) {
1303  this_point->neighbour_dist2 = numeric_limits<double>::max();
1304  // among all three shifts
1305  for (unsigned int ishift = 0; ishift < _nshift; ishift++) {
1306  circulator other = this_point->circ[ishift];
1307  // among points within CP_range
1308  for (unsigned i=0; i < CP_range; i++) {
1309  ++other;
1310  double dist2 = this_point->distance2(*other->point);
1311  if (dist2 < this_point->neighbour_dist2) {
1312  this_point->neighbour_dist2 = dist2;
1313  this_point->neighbour = other->point;
1314  }
1315  }
1316  }
1317  }
1318  _heap->update(_ID(this_point), this_point->neighbour_dist2);
1319  }
1320  this_point->review_flag = 0;
1321  }
1322 }
1323 unsigned int ClosestPair2D::insert(const Coord2D & new_coord) {
1324  assert(_available_points.size() > 0);
1325  Point * new_point = _available_points.top();
1326  _available_points.pop();
1327  new_point->coord = new_coord;
1328  _insert_into_search_tree(new_point);
1329  _deal_with_points_to_review();
1330  return _ID(new_point);
1331 }
1332 unsigned int ClosestPair2D::replace(unsigned int ID1, unsigned int ID2,
1333  const Coord2D & position) {
1334  Point * point_to_remove = & (_points[ID1]);
1335  _remove_from_search_tree(point_to_remove);
1336  point_to_remove = & (_points[ID2]);
1337  _remove_from_search_tree(point_to_remove);
1338  Point * new_point = _available_points.top();
1339  _available_points.pop();
1340  new_point->coord = position;
1341  _insert_into_search_tree(new_point);
1342  _deal_with_points_to_review();
1343  return _ID(new_point);
1344 }
1346  const std::vector<unsigned int> & IDs_to_remove,
1347  const std::vector<Coord2D> & new_positions,
1348  std::vector<unsigned int> & new_IDs) {
1349  for (unsigned int i = 0; i < IDs_to_remove.size(); i++) {
1350  _remove_from_search_tree(& (_points[IDs_to_remove[i]]));
1351  }
1352  new_IDs.resize(0);
1353  for (unsigned int i = 0; i < new_positions.size(); i++) {
1354  Point * new_point = _available_points.top();
1355  _available_points.pop();
1356  new_point->coord = new_positions[i];
1357  _insert_into_search_tree(new_point);
1358  new_IDs.push_back(_ID(new_point));
1359  }
1360  _deal_with_points_to_review();
1361 }
1363  _set_label(new_point, _review_heap_entry);
1364  new_point->neighbour_dist2 = numeric_limits<double>::max();
1365  unsigned int CP_range = min(_cp_search_range, size()-1);
1366  for (unsigned ishift = 0; ishift < _nshift; ishift++) {
1367  Shuffle new_shuffle;
1368  _point2shuffle(*new_point, new_shuffle, _shifts[ishift]);
1369  circulator new_circ = _trees[ishift]->insert(new_shuffle);
1370  new_point->circ[ishift] = new_circ;
1371  circulator right_edge = new_circ; right_edge++;
1372  circulator left_edge = new_circ;
1373  for (unsigned int i = 0; i < CP_range; i++) {left_edge--;}
1374  do {
1375  Point * left_point = left_edge->point;
1376  Point * right_point = right_edge->point;
1377  double new_dist2 = left_point->distance2(*new_point);
1378  if (new_dist2 < left_point->neighbour_dist2) {
1379  left_point->neighbour_dist2 = new_dist2;
1380  left_point->neighbour = new_point;
1381  _add_label(left_point, _review_heap_entry);
1382  }
1383  new_dist2 = new_point->distance2(*right_point);
1384  if (new_dist2 < new_point->neighbour_dist2) {
1385  new_point->neighbour_dist2 = new_dist2;
1386  new_point->neighbour = right_point;
1387  }
1388  if (left_point->neighbour == right_point) {
1389  _add_label(left_point, _review_neighbour);
1390  }
1391  right_edge++;
1392  } while (++left_edge != new_circ);
1393  }
1394 }
1396 #include<iostream>
1397 #include<sstream>
1398 #include<fstream>
1399 #include<cmath>
1400 #include<cstdlib>
1401 #include<cassert>
1402 #include<string>
1403 #include<set>
1404 FJCORE_BEGIN_NAMESPACE // defined in fastjet/internal/base.hh
1405 using namespace std;
1406 std::ostream * ClusterSequence::_fastjet_banner_ostr = &cout;
1408  if (_structure_shared_ptr){
1409  ClusterSequenceStructure* csi = dynamic_cast<ClusterSequenceStructure*>(_structure_shared_ptr.get());
1410  assert(csi != NULL);
1411  csi->set_associated_cs(NULL);
1412  if (_deletes_self_when_unused) {
1413  _structure_shared_ptr.set_count(_structure_shared_ptr.use_count()
1414  + _structure_use_count_after_construction);
1415  }
1416  }
1417 }
1419  assert(_deletes_self_when_unused);
1420  _deletes_self_when_unused = false;
1421 }
1423  const JetDefinition & jet_def_in,
1424  const bool & writeout_combinations) {
1425  _decant_options(jet_def_in, writeout_combinations);
1426  _initialise_and_run_no_decant();
1427 }
1429  _fill_initial_history();
1430  if (n_particles() == 0) return;
1431  if (_jet_algorithm == plugin_algorithm) {
1432  _plugin_activated = true;
1433  _jet_def.plugin()->run_clustering( (*this) );
1434  _plugin_activated = false;
1435  _update_structure_use_count();
1436  return;
1437  } else if (_jet_algorithm == ee_kt_algorithm ||
1438  _jet_algorithm == ee_genkt_algorithm) {
1439  _strategy = N2Plain;
1440  if (_jet_algorithm == ee_kt_algorithm) {
1441  assert(_Rparam > 2.0);
1442  _invR2 = 1.0;
1443  } else {
1444  if (_Rparam > pi) {
1445  // choose a value that ensures that back-to-back particles will
1446  // always recombine
1447  //_R2 = 4.0000000000001;
1448  _R2 = 2 * ( 3.0 + cos(_Rparam) );
1449  } else {
1450  _R2 = 2 * ( 1.0 - cos(_Rparam) );
1451  }
1452  _invR2 = 1.0/_R2;
1453  }
1454  _simple_N2_cluster_EEBriefJet();
1455  return;
1456  } else if (_jet_algorithm == undefined_jet_algorithm) {
1457  throw Error("A ClusterSequence cannot be created with an uninitialised JetDefinition");
1458  }
1459  if (_strategy == Best) {
1460  _strategy = _best_strategy();
1461 #ifdef __FJCORE_DROP_CGAL
1462  if (_strategy == NlnN) _strategy = N2MHTLazy25;
1463 #endif // __FJCORE_DROP_CGAL
1464  } else if (_strategy == BestFJ30) {
1465  int N = _jets.size();
1466  if (min(1.0,max(0.1,_Rparam)*3.3)*N <= 30) {
1467  _strategy = N2Plain;
1468  } else if (N > 6200/pow(_Rparam,2.0) && _jet_def.jet_algorithm() == cambridge_algorithm) {
1469  _strategy = NlnNCam;
1470 #ifndef __FJCORE_DROP_CGAL
1471  } else if ((N > 16000/pow(_Rparam,1.15) && _jet_def.jet_algorithm() != antikt_algorithm)
1472  || N > 35000/pow(_Rparam,1.15)) {
1473  _strategy = NlnN;
1474 #endif // __FJCORE_DROP_CGAL
1475  } else if (N <= 450) {
1476  _strategy = N2Tiled;
1477  } else {
1478  _strategy = N2MinHeapTiled;
1479  }
1480  }
1481  if (_Rparam >= twopi) {
1482  if ( _strategy == NlnN
1483  || _strategy == NlnN3pi
1484  || _strategy == NlnNCam
1485  || _strategy == NlnNCam2pi2R
1486  || _strategy == NlnNCam4pi) {
1487 #ifdef __FJCORE_DROP_CGAL
1488  _strategy = N2MinHeapTiled;
1489 #else
1490  _strategy = NlnN4pi;
1491 #endif
1492  }
1493  if (_jet_def.strategy() != Best && _strategy != _jet_def.strategy()) {
1494  ostringstream oss;
1495  oss << "Cluster strategy " << strategy_string(_jet_def.strategy())
1496  << " automatically changed to " << strategy_string()
1497  << " because the former is not supported for R = " << _Rparam
1498  << " >= 2pi";
1499  _changed_strategy_warning.warn(oss.str());
1500  }
1501  }
1502  if (_strategy == N2Plain) {
1503  this->_simple_N2_cluster_BriefJet();
1504  } else if (_strategy == N2Tiled) {
1505  this->_faster_tiled_N2_cluster();
1506  } else if (_strategy == N2MinHeapTiled) {
1507  this->_minheap_faster_tiled_N2_cluster();
1508  } else if (_strategy == N2MHTLazy9Alt) {
1509  _plugin_activated = true;
1510  LazyTiling9Alt tiling(*this);
1511  tiling.run();
1512  _plugin_activated = false;
1513  } else if (_strategy == N2MHTLazy25) {
1514  _plugin_activated = true;
1515  LazyTiling25 tiling(*this);
1516  tiling.run();
1517  _plugin_activated = false;
1518  } else if (_strategy == N2MHTLazy9) {
1519  _plugin_activated = true;
1520  LazyTiling9 tiling(*this);
1521  tiling.run();
1522  _plugin_activated = false;
1523  } else if (_strategy == N2MHTLazy9AntiKtSeparateGhosts) {
1524  throw Error("N2MHTLazy9AntiKtSeparateGhosts strategy not supported with FJCORE");
1525  } else if (_strategy == NlnN) {
1526  this->_delaunay_cluster();
1527  } else if (_strategy == NlnNCam) {
1528  this->_CP2DChan_cluster_2piMultD();
1529  } else if (_strategy == NlnN3pi || _strategy == NlnN4pi ) {
1530  this->_delaunay_cluster();
1531  } else if (_strategy == N3Dumb ) {
1532  this->_really_dumb_cluster();
1533  } else if (_strategy == N2PoorTiled) {
1534  this->_tiled_N2_cluster();
1535  } else if (_strategy == NlnNCam4pi) {
1536  this->_CP2DChan_cluster();
1537  } else if (_strategy == NlnNCam2pi2R) {
1538  this->_CP2DChan_cluster_2pi2R();
1539  } else {
1540  ostringstream err;
1541  err << "Unrecognised value for strategy: "<<_strategy;
1542  throw Error(err.str());
1543  }
1544 }
1545 bool ClusterSequence::_first_time = true;
1548  return "FastJet version "+string(fastjet_version)+" [fjcore]";
1549 }
1551  if (!_first_time) {return;}
1552  _first_time = false;
1553  ostream * ostr = _fastjet_banner_ostr;
1554  if (!ostr) return;
1555  (*ostr) << "#--------------------------------------------------------------------------\n";
1556  (*ostr) << "# FastJet release " << fastjet_version << " [fjcore]" << endl;
1557  (*ostr) << "# M. Cacciari, G.P. Salam and G. Soyez \n";
1558  (*ostr) << "# A software package for jet finding and analysis at colliders \n";
1559  (*ostr) << "# http://fastjet.fr \n";
1560  (*ostr) << "# \n";
1561  (*ostr) << "# Please cite EPJC72(2012)1896 [arXiv:1111.6097] if you use this package\n";
1562  (*ostr) << "# for scientific work and optionally PLB641(2006)57 [hep-ph/0512210]. \n";
1563  (*ostr) << "# \n";
1564  (*ostr) << "# FastJet is provided without warranty under the terms of the GNU GPLv2.\n";
1565  (*ostr) << "# It uses T. Chan's closest pair algorithm, S. Fortune's Voronoi code";
1566 #ifndef __FJCORE_DROP_CGAL
1567  (*ostr) << ",\n# CGAL ";
1568 #else
1569  (*ostr) << "\n# ";
1570 #endif // __FJCORE_DROP_CGAL
1571  (*ostr) << "and 3rd party plugin jet algorithms. See COPYING file for details.\n";
1572  (*ostr) << "#--------------------------------------------------------------------------\n";
1573  ostr->flush();
1574 }
1575 void ClusterSequence::_decant_options(const JetDefinition & jet_def_in,
1576  const bool & writeout_combinations) {
1577  _jet_def = jet_def_in;
1578  _writeout_combinations = writeout_combinations;
1579  _structure_shared_ptr.reset(new ClusterSequenceStructure(this));
1580  _decant_options_partial();
1581 }
1583  print_banner();
1584  _jet_algorithm = _jet_def.jet_algorithm();
1585  _Rparam = _jet_def.R(); _R2 = _Rparam*_Rparam; _invR2 = 1.0/_R2;
1586  _strategy = _jet_def.strategy();
1587  _plugin_activated = false;
1588  _update_structure_use_count(); // make sure it's correct already here
1589 }
1591  _jets.reserve(_jets.size()*2);
1592  _history.reserve(_jets.size()*2);
1593  _Qtot = 0;
1594  for (int i = 0; i < static_cast<int>(_jets.size()) ; i++) {
1596  element.parent1 = InexistentParent;
1597  element.parent2 = InexistentParent;
1598  element.child = Invalid;
1599  element.jetp_index = i;
1600  element.dij = 0.0;
1601  element.max_dij_so_far = 0.0;
1602  _history.push_back(element);
1603  _jet_def.recombiner()->preprocess(_jets[i]);
1604  _jets[i].set_cluster_hist_index(i);
1605  _set_structure_shared_ptr(_jets[i]);
1606  _Qtot += _jets[i].E();
1607  }
1608  _initial_n = _jets.size();
1609  _deletes_self_when_unused = false;
1610 }
1611 string ClusterSequence::strategy_string (Strategy strategy_in) const {
1612  string strategy;
1613  switch(strategy_in) {
1614  case NlnN:
1615  strategy = "NlnN"; break;
1616  case NlnN3pi:
1617  strategy = "NlnN3pi"; break;
1618  case NlnN4pi:
1619  strategy = "NlnN4pi"; break;
1620  case N2Plain:
1621  strategy = "N2Plain"; break;
1622  case N2Tiled:
1623  strategy = "N2Tiled"; break;
1624  case N2MinHeapTiled:
1625  strategy = "N2MinHeapTiled"; break;
1626  case N2PoorTiled:
1627  strategy = "N2PoorTiled"; break;
1628  case N2MHTLazy9:
1629  strategy = "N2MHTLazy9"; break;
1630  case N2MHTLazy9Alt:
1631  strategy = "N2MHTLazy9Alt"; break;
1632  case N2MHTLazy25:
1633  strategy = "N2MHTLazy25"; break;
1635  strategy = "N2MHTLazy9AntiKtSeparateGhosts"; break;
1636  case N3Dumb:
1637  strategy = "N3Dumb"; break;
1638  case NlnNCam4pi:
1639  strategy = "NlnNCam4pi"; break;
1640  case NlnNCam2pi2R:
1641  strategy = "NlnNCam2pi2R"; break;
1642  case NlnNCam:
1643  strategy = "NlnNCam"; break; // 2piMultD
1644  case plugin_strategy:
1645  strategy = "plugin strategy"; break;
1646  default:
1647  strategy = "Unrecognized";
1648  }
1649  return strategy;
1650 }
1652  const PseudoJet & jet) const {
1653  if (_jet_algorithm == kt_algorithm) {return jet.kt2();}
1654  else if (_jet_algorithm == cambridge_algorithm) {return 1.0;}
1655  else if (_jet_algorithm == antikt_algorithm) {
1656  double kt2=jet.kt2();
1657  return kt2 > 1e-300 ? 1.0/kt2 : 1e300;
1658  } else if (_jet_algorithm == genkt_algorithm) {
1659  double kt2 = jet.kt2();
1660  double p = jet_def().extra_param();
1661  if (p <= 0 && kt2 < 1e-300) kt2 = 1e-300; // dodgy safety check
1662  return pow(kt2, p);
1663  } else if (_jet_algorithm == cambridge_for_passive_algorithm) {
1664  double kt2 = jet.kt2();
1665  double lim = _jet_def.extra_param();
1666  if (kt2 < lim*lim && kt2 != 0.0) {
1667  return 1.0/kt2;
1668  } else {return 1.0;}
1669  } else {throw Error("Unrecognised jet algorithm");}
1670 }
1672  int N = _jets.size();
1673  double bounded_R = max(_Rparam, 0.1);
1674  if (N <= 30 || N <= 39.0/(bounded_R + 0.6)) {
1675  return N2Plain;
1676  }
1677  const static _Parabola N_Tiled_to_MHT_lowR (-45.4947,54.3528,44.6283);
1678  const static _Parabola L_MHT_to_MHTLazy9_lowR (0.677807,-1.05006,10.6994);
1679  const static _Parabola L_MHTLazy9_to_MHTLazy25_akt_lowR(0.169967,-0.512589,12.1572);
1680  const static _Parabola L_MHTLazy9_to_MHTLazy25_kt_lowR (0.16237,-0.484612,12.3373);
1681  const static _Parabola L_MHTLazy9_to_MHTLazy25_cam_lowR = L_MHTLazy9_to_MHTLazy25_kt_lowR;
1682  const static _Parabola L_MHTLazy25_to_NlnN_akt_lowR (0.0472051,-0.22043,15.9196);
1683  const static _Parabola L_MHTLazy25_to_NlnN_kt_lowR (0.118609,-0.326811,14.8287);
1684  const static _Parabola L_MHTLazy25_to_NlnN_cam_lowR (0.10119,-0.295748,14.3924);
1685  const static _Line L_Tiled_to_MHTLazy9_medR (-1.31304,7.29621);
1686  const static _Parabola L_MHTLazy9_to_MHTLazy25_akt_medR = L_MHTLazy9_to_MHTLazy25_akt_lowR;
1687  const static _Parabola L_MHTLazy9_to_MHTLazy25_kt_medR = L_MHTLazy9_to_MHTLazy25_kt_lowR;
1688  const static _Parabola L_MHTLazy9_to_MHTLazy25_cam_medR = L_MHTLazy9_to_MHTLazy25_cam_lowR;
1689  const static _Parabola L_MHTLazy25_to_NlnN_akt_medR = L_MHTLazy25_to_NlnN_akt_lowR;
1690  const static _Parabola L_MHTLazy25_to_NlnN_kt_medR = L_MHTLazy25_to_NlnN_kt_lowR;
1691  const static _Parabola L_MHTLazy25_to_NlnN_cam_medR = L_MHTLazy25_to_NlnN_cam_lowR;
1692  const static double N_Plain_to_MHTLazy9_largeR = 75;
1693  const static double N_MHTLazy9_to_MHTLazy25_akt_largeR = 700;
1694  const static double N_MHTLazy9_to_MHTLazy25_kt_largeR = 1000;
1695  const static double N_MHTLazy9_to_MHTLazy25_cam_largeR = 1000;
1696  const static double N_MHTLazy25_to_NlnN_akt_largeR = 100000;
1697  const static double N_MHTLazy25_to_NlnN_kt_largeR = 40000;
1698  const static double N_MHTLazy25_to_NlnN_cam_largeR = 15000;
1699  JetAlgorithm jet_algorithm;
1700  if (_jet_algorithm == genkt_algorithm) {
1701  double p = jet_def().extra_param();
1702  if (p < 0.0) jet_algorithm = antikt_algorithm;
1703  else jet_algorithm = kt_algorithm;
1704  } else if (_jet_algorithm == cambridge_for_passive_algorithm) {
1705  jet_algorithm = kt_algorithm;
1706  } else {
1707  jet_algorithm = _jet_algorithm;
1708  }
1709  if (bounded_R < 0.65) {
1710  if (N < N_Tiled_to_MHT_lowR(bounded_R)) return N2Tiled;
1711  double logN = log(double(N));
1712  if (logN < L_MHT_to_MHTLazy9_lowR(bounded_R)) return N2MinHeapTiled;
1713  else {
1714  if (jet_algorithm == antikt_algorithm){
1715  if (logN < L_MHTLazy9_to_MHTLazy25_akt_lowR(bounded_R)) return N2MHTLazy9;
1716  else if (logN < L_MHTLazy25_to_NlnN_akt_lowR(bounded_R)) return N2MHTLazy25;
1717  else return NlnN;
1718  } else if (jet_algorithm == kt_algorithm){
1719  if (logN < L_MHTLazy9_to_MHTLazy25_kt_lowR(bounded_R)) return N2MHTLazy9;
1720  else if (logN < L_MHTLazy25_to_NlnN_kt_lowR(bounded_R)) return N2MHTLazy25;
1721  else return NlnN;
1722  } else if (jet_algorithm == cambridge_algorithm) {
1723  if (logN < L_MHTLazy9_to_MHTLazy25_cam_lowR(bounded_R)) return N2MHTLazy9;
1724  else if (logN < L_MHTLazy25_to_NlnN_cam_lowR(bounded_R)) return N2MHTLazy25;
1725  else return NlnNCam;
1726  }
1727  }
1728  } else if (bounded_R < 0.5*pi) {
1729  double logN = log(double(N));
1730  if (logN < L_Tiled_to_MHTLazy9_medR(bounded_R)) return N2Tiled;
1731  else {
1732  if (jet_algorithm == antikt_algorithm){
1733  if (logN < L_MHTLazy9_to_MHTLazy25_akt_medR(bounded_R)) return N2MHTLazy9;
1734  else if (logN < L_MHTLazy25_to_NlnN_akt_medR(bounded_R)) return N2MHTLazy25;
1735  else return NlnN;
1736  } else if (jet_algorithm == kt_algorithm){
1737  if (logN < L_MHTLazy9_to_MHTLazy25_kt_medR(bounded_R)) return N2MHTLazy9;
1738  else if (logN < L_MHTLazy25_to_NlnN_kt_medR(bounded_R)) return N2MHTLazy25;
1739  else return NlnN;
1740  } else if (jet_algorithm == cambridge_algorithm) {
1741  if (logN < L_MHTLazy9_to_MHTLazy25_cam_medR(bounded_R)) return N2MHTLazy9;
1742  else if (logN < L_MHTLazy25_to_NlnN_cam_medR(bounded_R)) return N2MHTLazy25;
1743  else return NlnNCam;
1744  }
1745  }
1746  } else {
1747  if (N < N_Plain_to_MHTLazy9_largeR) return N2Plain;
1748  else {
1749  if (jet_algorithm == antikt_algorithm){
1750  if (N < N_MHTLazy9_to_MHTLazy25_akt_largeR) return N2MHTLazy9;
1751  else if (N < N_MHTLazy25_to_NlnN_akt_largeR) return N2MHTLazy25;
1752  else return NlnN;
1753  } else if (jet_algorithm == kt_algorithm){
1754  if (N < N_MHTLazy9_to_MHTLazy25_kt_largeR) return N2MHTLazy9;
1755  else if (N < N_MHTLazy25_to_NlnN_kt_largeR) return N2MHTLazy25;
1756  else return NlnN;
1757  } else if (jet_algorithm == cambridge_algorithm) {
1758  if (N < N_MHTLazy9_to_MHTLazy25_cam_largeR) return N2MHTLazy9;
1759  else if (N < N_MHTLazy25_to_NlnN_cam_largeR) return N2MHTLazy25;
1760  else return NlnNCam;
1761  }
1762  }
1763  }
1764  assert(0 && "Code should never reach here");
1765  return N2MHTLazy9;
1766 }
1767 ClusterSequence & ClusterSequence::operator=(const ClusterSequence & cs) {
1768  if (&cs != this) {
1769  _deletes_self_when_unused = false;
1770  transfer_from_sequence(cs);
1771  }
1772  return *this;
1773 }
1774 void ClusterSequence::transfer_from_sequence(const ClusterSequence & from_seq,
1775  const FunctionOfPseudoJet<PseudoJet> * action_on_jets){
1776  if (will_delete_self_when_unused())
1777  throw(Error("cannot use CS::transfer_from_sequence after a call to delete_self_when_unused()"));
1778  _jet_def = from_seq._jet_def ;
1779  _writeout_combinations = from_seq._writeout_combinations ;
1780  _initial_n = from_seq._initial_n ;
1781  _Rparam = from_seq._Rparam ;
1782  _R2 = from_seq._R2 ;
1783  _invR2 = from_seq._invR2 ;
1784  _strategy = from_seq._strategy ;
1785  _jet_algorithm = from_seq._jet_algorithm ;
1786  _plugin_activated = from_seq._plugin_activated ;
1787  if (action_on_jets)
1788  _jets = (*action_on_jets)(from_seq._jets);
1789  else
1790  _jets = from_seq._jets;
1791  _history = from_seq._history;
1792  _extras = from_seq._extras;
1793  if (_structure_shared_ptr) {
1794  if (_deletes_self_when_unused) throw Error("transfer_from_sequence cannot be used for a cluster sequence that deletes self when unused");
1795  ClusterSequenceStructure* csi = dynamic_cast<ClusterSequenceStructure*>(_structure_shared_ptr.get());
1796  assert(csi != NULL);
1797  csi->set_associated_cs(NULL);
1798  }
1799  _structure_shared_ptr.reset(new ClusterSequenceStructure(this));
1800  _update_structure_use_count();
1801  for (unsigned int i=0; i<_jets.size(); i++){
1802  _jets[i].set_cluster_hist_index(from_seq._jets[i].cluster_hist_index());
1803  _set_structure_shared_ptr(_jets[i]);
1804  }
1805 }
1807  int jet_i, int jet_j, double dij,
1808  const PseudoJet & newjet, int & newjet_k) {
1809  plugin_record_ij_recombination(jet_i, jet_j, dij, newjet_k);
1810  int tmp_index = _jets[newjet_k].cluster_hist_index();
1811  _jets[newjet_k] = newjet;
1812  _jets[newjet_k].set_cluster_hist_index(tmp_index);
1813  _set_structure_shared_ptr(_jets[newjet_k]);
1814 }
1815 vector<PseudoJet> ClusterSequence::inclusive_jets (const double ptmin) const{
1816  double dcut = ptmin*ptmin;
1817  int i = _history.size() - 1; // last jet
1818  vector<PseudoJet> jets_local;
1819  if (_jet_algorithm == kt_algorithm) {
1820  while (i >= 0) {
1821  if (_history[i].max_dij_so_far < dcut) {break;}
1822  if (_history[i].parent2 == BeamJet && _history[i].dij >= dcut) {
1823  // for beam jets
1824  int parent1 = _history[i].parent1;
1825  jets_local.push_back(_jets[_history[parent1].jetp_index]);}
1826  i--;
1827  }
1828  } else if (_jet_algorithm == cambridge_algorithm) {
1829  while (i >= 0) {
1830  if (_history[i].parent2 != BeamJet) {break;}
1831  int parent1 = _history[i].parent1;
1832  const PseudoJet & jet = _jets[_history[parent1].jetp_index];
1833  if (jet.perp2() >= dcut) {jets_local.push_back(jet);}
1834  i--;
1835  }
1836  } else if (_jet_algorithm == plugin_algorithm
1837  || _jet_algorithm == ee_kt_algorithm
1838  || _jet_algorithm == antikt_algorithm
1839  || _jet_algorithm == genkt_algorithm
1840  || _jet_algorithm == ee_genkt_algorithm
1841  || _jet_algorithm == cambridge_for_passive_algorithm) {
1842  while (i >= 0) {
1843  if (_history[i].parent2 == BeamJet) {
1844  int parent1 = _history[i].parent1;
1845  const PseudoJet & jet = _jets[_history[parent1].jetp_index];
1846  if (jet.perp2() >= dcut) {jets_local.push_back(jet);}
1847  }
1848  i--;
1849  }
1850  } else {throw Error("cs::inclusive_jets(...): Unrecognized jet algorithm");}
1851  return jets_local;
1852 }
1853 int ClusterSequence::n_exclusive_jets (const double dcut) const {
1854  int i = _history.size() - 1; // last jet
1855  while (i >= 0) {
1856  if (_history[i].max_dij_so_far <= dcut) {break;}
1857  i--;
1858  }
1859  int stop_point = i + 1;
1860  int njets = 2*_initial_n - stop_point;
1861  return njets;
1862 }
1863 vector<PseudoJet> ClusterSequence::exclusive_jets (const double dcut) const {
1864  int njets = n_exclusive_jets(dcut);
1865  return exclusive_jets(njets);
1866 }
1867 vector<PseudoJet> ClusterSequence::exclusive_jets (const int njets) const {
1868  if (njets > _initial_n) {
1869  ostringstream err;
1870  err << "Requested " << njets << " exclusive jets, but there were only "
1871  << _initial_n << " particles in the event";
1872  throw Error(err.str());
1873  }
1874  return exclusive_jets_up_to(njets);
1875 }
1876 vector<PseudoJet> ClusterSequence::exclusive_jets_up_to (const int njets) const {
1877  if (( _jet_def.jet_algorithm() != kt_algorithm) &&
1878  ( _jet_def.jet_algorithm() != cambridge_algorithm) &&
1879  ( _jet_def.jet_algorithm() != ee_kt_algorithm) &&
1880  (((_jet_def.jet_algorithm() != genkt_algorithm) &&
1881  (_jet_def.jet_algorithm() != ee_genkt_algorithm)) ||
1882  (_jet_def.extra_param() <0)) &&
1883  ((_jet_def.jet_algorithm() != plugin_algorithm) ||
1884  (!_jet_def.plugin()->exclusive_sequence_meaningful()))) {
1885  _exclusive_warnings.warn("dcut and exclusive jets for jet-finders other than kt, C/A or genkt with p>=0 should be interpreted with care.");
1886  }
1887  int stop_point = 2*_initial_n - njets;
1888  if (stop_point < _initial_n) stop_point = _initial_n;
1889  if (2*_initial_n != static_cast<int>(_history.size())) {
1890  ostringstream err;
1891  err << "2*_initial_n != _history.size() -- this endangers internal assumptions!\n";
1892  throw Error(err.str());
1893  }
1894  vector<PseudoJet> jets_local;
1895  for (unsigned int i = stop_point; i < _history.size(); i++) {
1896  int parent1 = _history[i].parent1;
1897  if (parent1 < stop_point) {
1898  jets_local.push_back(_jets[_history[parent1].jetp_index]);
1899  }
1900  int parent2 = _history[i].parent2;
1901  if (parent2 < stop_point && parent2 > 0) {
1902  jets_local.push_back(_jets[_history[parent2].jetp_index]);
1903  }
1904  }
1905  if (int(jets_local.size()) != min(_initial_n, njets)) {
1906  ostringstream err;
1907  err << "ClusterSequence::exclusive_jets: size of returned vector ("
1908  <<jets_local.size()<<") does not coincide with requested number of jets ("
1909  <<njets<<")";
1910  throw Error(err.str());
1911  }
1912  return jets_local;
1913 }
1914 double ClusterSequence::exclusive_dmerge (const int njets) const {
1915  assert(njets >= 0);
1916  if (njets >= _initial_n) {return 0.0;}
1917  return _history[2*_initial_n-njets-1].dij;
1918 }
1919 double ClusterSequence::exclusive_dmerge_max (const int njets) const {
1920  assert(njets >= 0);
1921  if (njets >= _initial_n) {return 0.0;}
1922  return _history[2*_initial_n-njets-1].max_dij_so_far;
1923 }
1924 std::vector<PseudoJet> ClusterSequence::exclusive_subjets
1925  (const PseudoJet & jet, const double dcut) const {
1926  set<const history_element*> subhist;
1927  get_subhist_set(subhist, jet, dcut, 0);
1928  vector<PseudoJet> subjets;
1929  subjets.reserve(subhist.size());
1930  for (set<const history_element*>::iterator elem = subhist.begin();
1931  elem != subhist.end(); elem++) {
1932  subjets.push_back(_jets[(*elem)->jetp_index]);
1933  }
1934  return subjets;
1935 }
1936 int ClusterSequence::n_exclusive_subjets(const PseudoJet & jet,
1937  const double dcut) const {
1938  set<const history_element*> subhist;
1939  get_subhist_set(subhist, jet, dcut, 0);
1940  return subhist.size();
1941 }
1942 std::vector<PseudoJet> ClusterSequence::exclusive_subjets
1943  (const PseudoJet & jet, int nsub) const {
1944  vector<PseudoJet> subjets = exclusive_subjets_up_to(jet, nsub);
1945  if (int(subjets.size()) < nsub) {
1946  ostringstream err;
1947  err << "Requested " << nsub << " exclusive subjets, but there were only "
1948  << subjets.size() << " particles in the jet";
1949  throw Error(err.str());
1950  }
1951  return subjets;
1952 }
1953 std::vector<PseudoJet> ClusterSequence::exclusive_subjets_up_to
1954  (const PseudoJet & jet, int nsub) const {
1955  set<const history_element*> subhist;
1956  vector<PseudoJet> subjets;
1957  if (nsub < 0) throw Error("Requested a negative number of subjets. This is nonsensical.");
1958  if (nsub == 0) return subjets;
1959  get_subhist_set(subhist, jet, -1.0, nsub);
1960  subjets.reserve(subhist.size());
1961  for (set<const history_element*>::iterator elem = subhist.begin();
1962  elem != subhist.end(); elem++) {
1963  subjets.push_back(_jets[(*elem)->jetp_index]);
1964  }
1965  return subjets;
1966 }
1967 double ClusterSequence::exclusive_subdmerge(const PseudoJet & jet, int nsub) const {
1968  set<const history_element*> subhist;
1969  get_subhist_set(subhist, jet, -1.0, nsub);
1970  set<const history_element*>::iterator highest = subhist.end();
1971  highest--;
1972  return (*highest)->dij;
1973 }
1974 double ClusterSequence::exclusive_subdmerge_max(const PseudoJet & jet, int nsub) const {
1975  set<const history_element*> subhist;
1976  get_subhist_set(subhist, jet, -1.0, nsub);
1977  set<const history_element*>::iterator highest = subhist.end();
1978  highest--;
1979  return (*highest)->max_dij_so_far;
1980 }
1981 void ClusterSequence::get_subhist_set(set<const history_element*> & subhist,
1982  const PseudoJet & jet,
1983  double dcut, int maxjet) const {
1984  assert(contains(jet));
1985  subhist.clear();
1986  subhist.insert(&(_history[jet.cluster_hist_index()]));
1987  int njet = 1;
1988  while (true) {
1989  set<const history_element*>::iterator highest = subhist.end();
1990  assert (highest != subhist.begin());
1991  highest--;
1992  const history_element* elem = *highest;
1993  if (njet == maxjet) break;
1994  if (elem->parent1 < 0) break;
1995  if (elem->max_dij_so_far <= dcut) break;
1996  subhist.erase(highest);
1997  subhist.insert(&(_history[elem->parent1]));
1998  subhist.insert(&(_history[elem->parent2]));
1999  njet++;
2000  }
2001 }
2002 bool ClusterSequence::object_in_jet(const PseudoJet & object,
2003  const PseudoJet & jet) const {
2004  assert(contains(object) && contains(jet));
2005  const PseudoJet * this_object = &object;
2006  const PseudoJet * childp;
2007  while(true) {
2008  if (this_object->cluster_hist_index() == jet.cluster_hist_index()) {
2009  return true;
2010  } else if (has_child(*this_object, childp)) {
2011  this_object = childp;
2012  } else {
2013  return false;
2014  }
2015  }
2016 }
2017 bool ClusterSequence::has_parents(const PseudoJet & jet, PseudoJet & parent1,
2018  PseudoJet & parent2) const {
2019  const history_element & hist = _history[jet.cluster_hist_index()];
2020  assert ((hist.parent1 >= 0 && hist.parent2 >= 0) ||
2021  (hist.parent1 < 0 && hist.parent2 < 0));
2022  if (hist.parent1 < 0) {
2023  parent1 = PseudoJet(0.0,0.0,0.0,0.0);
2024  parent2 = parent1;
2025  return false;
2026  } else {
2027  parent1 = _jets[_history[hist.parent1].jetp_index];
2028  parent2 = _jets[_history[hist.parent2].jetp_index];
2029  if (parent1.perp2() < parent2.perp2()) std::swap(parent1,parent2);
2030  return true;
2031  }
2032 }
2033 bool ClusterSequence::has_child(const PseudoJet & jet, PseudoJet & child) const {
2034  const PseudoJet * childp;
2035  bool res = has_child(jet, childp);
2036  if (res) {
2037  child = *childp;
2038  return true;
2039  } else {
2040  child = PseudoJet(0.0,0.0,0.0,0.0);
2041  return false;
2042  }
2043 }
2044 bool ClusterSequence::has_child(const PseudoJet & jet, const PseudoJet * & childp) const {
2045  const history_element & hist = _history[jet.cluster_hist_index()];
2046  if (hist.child >= 0 && _history[hist.child].jetp_index >= 0) {
2047  childp = &(_jets[_history[hist.child].jetp_index]);
2048  return true;
2049  } else {
2050  childp = NULL;
2051  return false;
2052  }
2053 }
2054 bool ClusterSequence::has_partner(const PseudoJet & jet,
2055  PseudoJet & partner) const {
2056  const history_element & hist = _history[jet.cluster_hist_index()];
2057  if (hist.child >= 0 && _history[hist.child].parent2 >= 0) {
2058  const history_element & child_hist = _history[hist.child];
2059  if (child_hist.parent1 == jet.cluster_hist_index()) {
2060  partner = _jets[_history[child_hist.parent2].jetp_index];
2061  } else {
2062  partner = _jets[_history[child_hist.parent1].jetp_index];
2063  }
2064  return true;
2065  } else {
2066  partner = PseudoJet(0.0,0.0,0.0,0.0);
2067  return false;
2068  }
2069 }
2070 vector<PseudoJet> ClusterSequence::constituents (const PseudoJet & jet) const {
2071  vector<PseudoJet> subjets;
2072  add_constituents(jet, subjets);
2073  return subjets;
2074 }
2075 void ClusterSequence::print_jets_for_root(const std::vector<PseudoJet> & jets_in,
2076  ostream & ostr) const {
2077  for (unsigned i = 0; i < jets_in.size(); i++) {
2078  ostr << i << " "
2079  << jets_in[i].px() << " "
2080  << jets_in[i].py() << " "
2081  << jets_in[i].pz() << " "
2082  << jets_in[i].E() << endl;
2083  vector<PseudoJet> cst = constituents(jets_in[i]);
2084  for (unsigned j = 0; j < cst.size() ; j++) {
2085  ostr << " " << j << " "
2086  << cst[j].rap() << " "
2087  << cst[j].phi() << " "
2088  << cst[j].perp() << endl;
2089  }
2090  ostr << "#END" << endl;
2091  }
2092 }
2093 void ClusterSequence::print_jets_for_root(const std::vector<PseudoJet> & jets_in,
2094  const std::string & filename,
2095  const std::string & comment ) const {
2096  std::ofstream ostr(filename.c_str());
2097  if (comment != "") ostr << "# " << comment << endl;
2098  print_jets_for_root(jets_in, ostr);
2099 }
2101  const vector<PseudoJet> & jets_in) const {
2102  vector<int> indices(n_particles());
2103  for (unsigned ipart = 0; ipart < n_particles(); ipart++)
2104  indices[ipart] = -1;
2105  for (unsigned ijet = 0; ijet < jets_in.size(); ijet++) {
2106  vector<PseudoJet> jet_constituents(constituents(jets_in[ijet]));
2107  for (unsigned ip = 0; ip < jet_constituents.size(); ip++) {
2108  unsigned iclust = jet_constituents[ip].cluster_hist_index();
2109  unsigned ipart = history()[iclust].jetp_index;
2110  indices[ipart] = ijet;
2111  }
2112  }
2113  return indices;
2114 }
2116  const PseudoJet & jet, vector<PseudoJet> & subjet_vector) const {
2117  int i = jet.cluster_hist_index();
2118  int parent1 = _history[i].parent1;
2119  int parent2 = _history[i].parent2;
2120  if (parent1 == InexistentParent) {
2121  subjet_vector.push_back(_jets[i]);
2122  return;
2123  }
2124  add_constituents(_jets[_history[parent1].jetp_index], subjet_vector);
2125  if (parent2 != BeamJet) {
2126  add_constituents(_jets[_history[parent2].jetp_index], subjet_vector);
2127  }
2128 }
2130  const int parent1,
2131  const int parent2, const int jetp_index,
2132  const double dij) {
2134  element.parent1 = parent1;
2135  element.parent2 = parent2;
2136  element.jetp_index = jetp_index;
2137  element.child = Invalid;
2138  element.dij = dij;
2139  element.max_dij_so_far = max(dij,_history[_history.size()-1].max_dij_so_far);
2140  _history.push_back(element);
2141  int local_step = _history.size()-1;
2142  assert(parent1 >= 0);
2143  if (_history[parent1].child != Invalid){
2144  throw InternalError("trying to recomine an object that has previsously been recombined");
2145  }
2146  _history[parent1].child = local_step;
2147  if (parent2 >= 0) {
2148  if (_history[parent2].child != Invalid){
2149  throw InternalError("trying to recomine an object that has previsously been recombined");
2150  }
2151  _history[parent2].child = local_step;
2152  }
2153  if (jetp_index != Invalid) {
2154  assert(jetp_index >= 0);
2155  _jets[jetp_index].set_cluster_hist_index(local_step);
2156  _set_structure_shared_ptr(_jets[jetp_index]);
2157  }
2158  if (_writeout_combinations) {
2159  cout << local_step << ": "
2160  << parent1 << " with " << parent2
2161  << "; y = "<< dij<<endl;
2162  }
2163 }
2165  valarray<int> lowest_constituent(_history.size());
2166  int hist_n = _history.size();
2167  lowest_constituent = hist_n; // give it a large number
2168  for (int i = 0; i < hist_n; i++) {
2169  lowest_constituent[i] = min(lowest_constituent[i],i);
2170  if (_history[i].child > 0) lowest_constituent[_history[i].child]
2171  = min(lowest_constituent[_history[i].child],lowest_constituent[i]);
2172  }
2173  valarray<bool> extracted(_history.size()); extracted = false;
2174  vector<int> unique_tree;
2175  unique_tree.reserve(_history.size());
2176  for (unsigned i = 0; i < n_particles(); i++) {
2177  if (!extracted[i]) {
2178  unique_tree.push_back(i);
2179  extracted[i] = true;
2180  _extract_tree_children(i, extracted, lowest_constituent, unique_tree);
2181  }
2182  }
2183  return unique_tree;
2184 }
2186  int position,
2187  valarray<bool> & extracted,
2188  const valarray<int> & lowest_constituent,
2189  vector<int> & unique_tree) const {
2190  if (!extracted[position]) {
2191  _extract_tree_parents(position,extracted,lowest_constituent,unique_tree);
2192  }
2193  int child = _history[position].child;
2194  if (child >= 0) _extract_tree_children(child,extracted,lowest_constituent,unique_tree);
2195 }
2196 vector<PseudoJet> ClusterSequence::unclustered_particles() const {
2197  vector<PseudoJet> unclustered;
2198  for (unsigned i = 0; i < n_particles() ; i++) {
2199  if (_history[i].child == Invalid)
2200  unclustered.push_back(_jets[_history[i].jetp_index]);
2201  }
2202  return unclustered;
2203 }
2204 vector<PseudoJet> ClusterSequence::childless_pseudojets() const {
2205  vector<PseudoJet> unclustered;
2206  for (unsigned i = 0; i < _history.size() ; i++) {
2207  if ((_history[i].child == Invalid) && (_history[i].parent2 != BeamJet))
2208  unclustered.push_back(_jets[_history[i].jetp_index]);
2209  }
2210  return unclustered;
2211 }
2212 bool ClusterSequence::contains(const PseudoJet & jet) const {
2213  return jet.cluster_hist_index() >= 0
2214  && jet.cluster_hist_index() < int(_history.size())
2215  && jet.has_valid_cluster_sequence()
2216  && jet.associated_cluster_sequence() == this;
2217 }
2219  int position,
2220  valarray<bool> & extracted,
2221  const valarray<int> & lowest_constituent,
2222  vector<int> & unique_tree) const {
2223  if (!extracted[position]) {
2224  int parent1 = _history[position].parent1;
2225  int parent2 = _history[position].parent2;
2226  if (parent1 >= 0 && parent2 >= 0) {
2227  if (lowest_constituent[parent1] > lowest_constituent[parent2])
2228  std::swap(parent1, parent2);
2229  }
2230  if (parent1 >= 0 && !extracted[parent1])
2231  _extract_tree_parents(parent1,extracted,lowest_constituent,unique_tree);
2232  if (parent2 >= 0 && !extracted[parent2])
2233  _extract_tree_parents(parent2,extracted,lowest_constituent,unique_tree);
2234  unique_tree.push_back(position);
2235  extracted[position] = true;
2236  }
2237 }
2239  const int jet_i, const int jet_j,
2240  const double dij,
2241  int & newjet_k) {
2242  PseudoJet newjet(false);
2243  _jet_def.recombiner()->recombine(_jets[jet_i], _jets[jet_j], newjet);
2244  _jets.push_back(newjet);
2245  newjet_k = _jets.size()-1;
2246  int newstep_k = _history.size();
2247  _jets[newjet_k].set_cluster_hist_index(newstep_k);
2248  int hist_i = _jets[jet_i].cluster_hist_index();
2249  int hist_j = _jets[jet_j].cluster_hist_index();
2250  _add_step_to_history(min(hist_i, hist_j), max(hist_i,hist_j),
2251  newjet_k, dij);
2252 }
2254  const int jet_i, const double diB) {
2255  _add_step_to_history(_jets[jet_i].cluster_hist_index(),BeamJet,
2256  Invalid, diB);
2257 }
2260  j.set_structure_shared_ptr(_structure_shared_ptr);
2261  _update_structure_use_count();
2262 }
2264  _structure_use_count_after_construction = _structure_shared_ptr.use_count();
2265 }
2267  int new_count = _structure_shared_ptr.use_count() - _structure_use_count_after_construction;
2268  if (new_count <= 0) {
2269  throw Error("delete_self_when_unused may only be called if at least one object outside the CS (e.g. a jet) is already associated with the CS");
2270  }
2271  _structure_shared_ptr.set_count(new_count);
2272  _deletes_self_when_unused = true;
2273 }
2275 #include<limits>
2276 #include<vector>
2277 #include<cmath>
2278 #include<iostream>
2279 FJCORE_BEGIN_NAMESPACE // defined in fastjet/internal/base.hh
2280 using namespace std;
2281 namespace Private {
2282  class MirrorInfo{
2283  public:
2284  int orig, mirror;
2285  MirrorInfo(int a, int b) : orig(a), mirror(b) {}
2286  MirrorInfo() : orig(0), mirror(0) {} // set dummy values to keep static code checkers happy
2287  };
2288  bool make_mirror(Coord2D & point, double Dlim) {
2289  if (point.y < Dlim) {point.y += twopi; return true;}
2290  if (twopi-point.y < Dlim) {point.y -= twopi; return true;}
2291  return false;
2292  }
2293 }
2294 using namespace Private;
2296  unsigned int n = _initial_n;
2297  vector<MirrorInfo> coordIDs(2*n); // coord IDs of a given jetID
2298  vector<int> jetIDs(2*n); // jet ID for a given coord ID
2299  vector<Coord2D> coords(2*n); // our coordinates (and copies)
2300  double Dlim4mirror = min(Dlim,pi);
2301  double minrap = numeric_limits<double>::max();
2302  double maxrap = -minrap;
2303  int coord_index = -1;
2304  int n_active = 0;
2305  for (unsigned jet_i = 0; jet_i < _jets.size(); jet_i++) {
2306  if (_history[_jets[jet_i].cluster_hist_index()].child != Invalid ||
2307  (_jets[jet_i].E() == abs(_jets[jet_i].pz()) &&
2308  _jets[jet_i].perp2() == 0.0)
2309  ) {continue;}
2310  n_active++;
2311  coordIDs[jet_i].orig = ++coord_index;
2312  coords[coord_index] = Coord2D(_jets[jet_i].rap(), _jets[jet_i].phi_02pi());
2313  jetIDs[coord_index] = jet_i;
2314  minrap = min(coords[coord_index].x,minrap);
2315  maxrap = max(coords[coord_index].x,maxrap);
2316  Coord2D mirror_point(coords[coord_index]);
2317  if (make_mirror(mirror_point, Dlim4mirror)) {
2318  coordIDs[jet_i].mirror = ++coord_index;
2319  coords[coord_index] = mirror_point;
2320  jetIDs[coord_index] = jet_i;
2321  } else {
2322  coordIDs[jet_i].mirror = Invalid;
2323  }
2324  }
2325  coords.resize(coord_index+1);
2326  Coord2D left_edge(minrap-1.0, -3.15); // a security margin below -pi
2327  Coord2D right_edge(maxrap+1.0, 9.45); // a security margin above 3*pi
2328  ClosestPair2D cp(coords, left_edge, right_edge);
2329  vector<Coord2D> new_points(2);
2330  vector<unsigned int> cIDs_to_remove(4);
2331  vector<unsigned int> new_cIDs(2);
2332  do {
2333  unsigned int cID1, cID2;
2334  double distance2;
2335  cp.closest_pair(cID1,cID2,distance2);
2336  if (distance2 > Dlim*Dlim) {break;}
2337  distance2 *= _invR2;
2338  int jet_i = jetIDs[cID1];
2339  int jet_j = jetIDs[cID2];
2340  assert (jet_i != jet_j); // to catch issue of recombining with mirror point
2341  int newjet_k;
2342  _do_ij_recombination_step(jet_i, jet_j, distance2, newjet_k);
2343  if (--n_active == 1) {break;}
2344  cIDs_to_remove.resize(0);
2345  cIDs_to_remove.push_back(coordIDs[jet_i].orig);
2346  cIDs_to_remove.push_back(coordIDs[jet_j].orig);
2347  if (coordIDs[jet_i].mirror != Invalid)
2348  cIDs_to_remove.push_back(coordIDs[jet_i].mirror);
2349  if (coordIDs[jet_j].mirror != Invalid)
2350  cIDs_to_remove.push_back(coordIDs[jet_j].mirror);
2351  Coord2D new_point(_jets[newjet_k].rap(),_jets[newjet_k].phi_02pi());
2352  new_points.resize(0);
2353  new_points.push_back(new_point);
2354  if (make_mirror(new_point, Dlim4mirror)) new_points.push_back(new_point); //< same warning as before concerning the mirroring
2355  cp.replace_many(cIDs_to_remove, new_points, new_cIDs);
2356  coordIDs[newjet_k].orig = new_cIDs[0];
2357  jetIDs[new_cIDs[0]] = newjet_k;
2358  if (new_cIDs.size() == 2) {
2359  coordIDs[newjet_k].mirror = new_cIDs[1];
2360  jetIDs[new_cIDs[1]] = newjet_k;
2361  } else {coordIDs[newjet_k].mirror = Invalid;}
2362  } while(true);
2363 }
2365  if (_jet_algorithm != cambridge_algorithm) throw Error("CP2DChan clustering method called for a jet-finder that is not the cambridge algorithm");
2366  _CP2DChan_limited_cluster(_Rparam);
2367  _do_Cambridge_inclusive_jets();
2368 }
2370  if (_Rparam >= 0.39) {
2371  _CP2DChan_limited_cluster(min(_Rparam/2,0.3));
2372  }
2373  _CP2DChan_cluster_2pi2R ();
2374 }
2376  if (_jet_algorithm != cambridge_algorithm) throw Error("_CP2DChan_cluster called for a jet-finder that is not the cambridge algorithm");
2377  unsigned int n = _jets.size();
2378  vector<MirrorInfo> coordIDs(2*n); // link from original to mirror indices
2379  vector<int> jetIDs(2*n); // link from mirror to original indices
2380  vector<Coord2D> coords(2*n); // our coordinates (and copies)
2381  double minrap = numeric_limits<double>::max();
2382  double maxrap = -minrap;
2383  int coord_index = 0;
2384  for (unsigned i = 0; i < n; i++) {
2385  if (_jets[i].E() == abs(_jets[i].pz()) && _jets[i].perp2() == 0.0) {
2386  coordIDs[i] = MirrorInfo(BeamJet,BeamJet);
2387  } else {
2388  coordIDs[i].orig = coord_index;
2389  coordIDs[i].mirror = coord_index+1;
2390  coords[coord_index] = Coord2D(_jets[i].rap(), _jets[i].phi_02pi());
2391  coords[coord_index+1] = Coord2D(_jets[i].rap(), _jets[i].phi_02pi()+twopi);
2392  jetIDs[coord_index] = i;
2393  jetIDs[coord_index+1] = i;
2394  minrap = min(coords[coord_index].x,minrap);
2395  maxrap = max(coords[coord_index].x,maxrap);
2396  coord_index += 2;
2397  }
2398  }
2399  for (unsigned i = n; i < 2*n; i++) {coordIDs[i].orig = Invalid;}
2400  coords.resize(coord_index);
2401  Coord2D left_edge(minrap-1.0, 0.0);
2402  Coord2D right_edge(maxrap+1.0, 2*twopi);
2403  ClosestPair2D cp(coords, left_edge, right_edge);
2404  vector<Coord2D> new_points(2);
2405  vector<unsigned int> cIDs_to_remove(4);
2406  vector<unsigned int> new_cIDs(2);
2407  do {
2408  unsigned int cID1, cID2;
2409  double distance2;
2410  cp.closest_pair(cID1,cID2,distance2);
2411  distance2 *= _invR2;
2412  if (distance2 > 1.0) {break;}
2413  int jet_i = jetIDs[cID1];
2414  int jet_j = jetIDs[cID2];
2415  assert (jet_i != jet_j); // to catch issue of recombining with mirror point
2416  int newjet_k;
2417  _do_ij_recombination_step(jet_i, jet_j, distance2, newjet_k);
2418  cIDs_to_remove[0] = coordIDs[jet_i].orig;
2419  cIDs_to_remove[1] = coordIDs[jet_i].mirror;
2420  cIDs_to_remove[2] = coordIDs[jet_j].orig;
2421  cIDs_to_remove[3] = coordIDs[jet_j].mirror;
2422  new_points[0] = Coord2D(_jets[newjet_k].rap(),_jets[newjet_k].phi_02pi());
2423  new_points[1] = Coord2D(_jets[newjet_k].rap(),_jets[newjet_k].phi_02pi()+twopi);
2424  new_cIDs[0] = cp.replace(cIDs_to_remove[0], cIDs_to_remove[2], new_points[0]);
2425  new_cIDs[1] = cp.replace(cIDs_to_remove[1], cIDs_to_remove[3], new_points[1]);
2426  coordIDs[jet_i].orig = Invalid;
2427  coordIDs[jet_j].orig = Invalid;
2428  coordIDs[newjet_k] = MirrorInfo(new_cIDs[0], new_cIDs[1]);
2429  jetIDs[new_cIDs[0]] = newjet_k;
2430  jetIDs[new_cIDs[1]] = newjet_k;
2431  n--;
2432  if (n == 1) {break;}
2433  } while(true);
2434  _do_Cambridge_inclusive_jets();
2435 }
2437  unsigned int n = _history.size();
2438  for (unsigned int hist_i = 0; hist_i < n; hist_i++) {
2439  if (_history[hist_i].child == Invalid) {
2440  _do_iB_recombination_step(_history[hist_i].jetp_index, 1.0);
2441  }
2442  }
2443 }
2445 #include<iostream>
2446 #include<sstream>
2447 #include<cmath>
2448 #include <cstdlib>
2449 #include<cassert>
2450 #include<memory>
2451 #ifndef __FJCORE_DROP_CGAL // in case we do not have the code for CGAL
2452 #include "fastjet/internal/Dnn4piCylinder.hh"
2453 #include "fastjet/internal/Dnn3piCylinder.hh"
2454 #include "fastjet/internal/Dnn2piCylinder.hh"
2455 #endif // __FJCORE_DROP_CGAL
2456 FJCORE_BEGIN_NAMESPACE // defined in fastjet/internal/base.hh
2457 using namespace std;
2459  int n = _jets.size();
2460  vector<EtaPhi> points(n); // recall EtaPhi is just a typedef'd pair<double>
2461  for (int i = 0; i < n; i++) {
2462  points[i] = EtaPhi(_jets[i].rap(),_jets[i].phi_02pi());
2463  points[i].sanitize(); // make sure things are in the right range
2464  }
2466  const bool verbose = false;
2467 #ifndef __FJCORE_DROP_CGAL // strategy = NlnN* are not supported if we drop CGAL...
2468  bool ignore_nearest_is_mirror = (_Rparam < twopi);
2469  if (_strategy == NlnN4pi) {
2470  DNN.reset(new Dnn4piCylinder(points,verbose));
2471  } else if (_strategy == NlnN3pi) {
2472  DNN.reset(new Dnn3piCylinder(points,ignore_nearest_is_mirror,verbose));
2473  } else if (_strategy == NlnN) {
2474  DNN.reset(new Dnn2piCylinder(points,ignore_nearest_is_mirror,verbose));
2475  } else
2476 #else
2477  if (_strategy == NlnN4pi || _strategy == NlnN3pi || _strategy == NlnN) {
2478  ostringstream err;
2479  err << "ERROR: Requested strategy "<<strategy_string()<<" but it is not"<<endl;
2480  err << " supported because FastJet was compiled without CGAL"<<endl;
2481  throw Error(err.str());
2482  } else
2483 #endif // __FJCORE_DROP_CGAL
2484  {
2485  assert(false);
2486  }
2487  DistMap DijMap;
2488  for (int ii = 0; ii < n; ii++) {
2489  _add_ktdistance_to_map(ii, DijMap, DNN.get());
2490  }
2491  for (int i=0;i<n;i++) {
2492  TwoVertices SmallestDijPair;
2493  int jet_i, jet_j;
2494  double SmallestDij;
2495  bool Valid2;
2496  bool recombine_with_beam;
2497  do {
2498  SmallestDij = DijMap.begin()->first;
2499  SmallestDijPair = DijMap.begin()->second;
2500  jet_i = SmallestDijPair.first;
2501  jet_j = SmallestDijPair.second;
2502  if (verbose) cout << "CS_Delaunay found recombination candidate: " << jet_i << " " << jet_j << " " << SmallestDij << endl; // GPS debugging
2503  DijMap.erase(DijMap.begin());
2504  recombine_with_beam = (jet_j == BeamJet);
2505  if (!recombine_with_beam) {Valid2 = DNN->Valid(jet_j);}
2506  else {Valid2 = true;}
2507  if (verbose) cout << "CS_Delaunay validities i & j: " << DNN->Valid(jet_i) << " " << Valid2 << endl;
2508  } while ( !DNN->Valid(jet_i) || !Valid2);
2509  if (! recombine_with_beam) {
2510  int nn; // will be index of new jet
2511  if (verbose) cout << "CS_Delaunay call _do_ij_recomb: " << jet_i << " " << jet_j << " " << SmallestDij << endl; // GPS debug
2512  _do_ij_recombination_step(jet_i, jet_j, SmallestDij, nn);
2513  EtaPhi newpoint(_jets[nn].rap(), _jets[nn].phi_02pi());
2514  newpoint.sanitize(); // make sure it is in correct range
2515  points.push_back(newpoint);
2516  } else {
2517  if (verbose) cout << "CS_Delaunay call _do_iB_recomb: " << jet_i << " " << SmallestDij << endl; // GPS debug
2518  _do_iB_recombination_step(jet_i, SmallestDij);
2519  }
2520  if (i == n-1) {break;}
2521  vector<int> updated_neighbours;
2522  if (! recombine_with_beam) {
2523  int point3;
2524  DNN->RemoveCombinedAddCombination(jet_i, jet_j,
2525  points[points.size()-1], point3,
2526  updated_neighbours);
2527  if (static_cast<unsigned int> (point3) != points.size()-1) {
2528  throw Error("INTERNAL ERROR: point3 != points.size()-1");}
2529  } else {
2530  DNN->RemovePoint(jet_i, updated_neighbours);
2531  }
2532  vector<int>::iterator it = updated_neighbours.begin();
2533  for (; it != updated_neighbours.end(); ++it) {
2534  int ii = *it;
2535  _add_ktdistance_to_map(ii, DijMap, DNN.get());
2536  }
2537  } // end clustering loop
2538 }
2540  const int ii,
2541  DistMap & DijMap,
2542  const DynamicNearestNeighbours * DNN) {
2543  double yiB = jet_scale_for_algorithm(_jets[ii]);
2544  if (yiB == 0.0) {
2545  DijMap.insert(DijEntry(yiB, TwoVertices(ii,-1)));
2546  } else {
2547  double DeltaR2 = DNN->NearestNeighbourDistance(ii) * _invR2;
2548  if (DeltaR2 > 1.0) {
2549  DijMap.insert(DijEntry(yiB, TwoVertices(ii,-1)));
2550  } else {
2551  double kt2i = jet_scale_for_algorithm(_jets[ii]);
2552  int jj = DNN->NearestNeighbourIndex(ii);
2553  if (kt2i <= jet_scale_for_algorithm(_jets[jj])) {
2554  double dij = DeltaR2 * kt2i;
2555  DijMap.insert(DijEntry(dij, TwoVertices(ii,jj)));
2556  }
2557  }
2558  }
2559 }
2561 #include<iostream>
2562 #include<cmath>
2563 #include <cstdlib>
2564 #include<cassert>
2565 FJCORE_BEGIN_NAMESPACE // defined in fastjet/internal/base.hh
2566 using namespace std;
2568  vector<PseudoJet *> jetsp(_jets.size());
2569  vector<int> indices(_jets.size());
2570  for (size_t i = 0; i<_jets.size(); i++) {
2571  jetsp[i] = & _jets[i];
2572  indices[i] = i;
2573  }
2574  for (int n = jetsp.size(); n > 0; n--) {
2575  int ii, jj;
2576  double ymin = jet_scale_for_algorithm(*(jetsp[0]));
2577  ii = 0; jj = -2;
2578  for (int i = 0; i < n; i++) {
2579  double yiB = jet_scale_for_algorithm(*(jetsp[i]));
2580  if (yiB < ymin) {
2581  ymin = yiB; ii = i; jj = -2;}
2582  }
2583  for (int i = 0; i < n-1; i++) {
2584  for (int j = i+1; j < n; j++) {
2585  //double y = jetsp[i]->kt_distance(*jetsp[j])*_invR2;
2586  double y = min(jet_scale_for_algorithm(*(jetsp[i])),
2587  jet_scale_for_algorithm(*(jetsp[j])))
2588  * jetsp[i]->plain_distance(*jetsp[j])*_invR2;
2589  if (y < ymin) {ymin = y; ii = i; jj = j;}
2590  }
2591  }
2592  int newn = 2*jetsp.size() - n;
2593  if (jj >= 0) {
2594  int nn; // new jet index
2595  _do_ij_recombination_step(jetsp[ii]-&_jets[0],
2596  jetsp[jj]-&_jets[0], ymin, nn);
2597  jetsp[ii] = &_jets[nn];
2598  jetsp[jj] = jetsp[n-1];
2599  indices[ii] = newn;
2600  indices[jj] = indices[n-1];
2601  } else {
2602  _do_iB_recombination_step(jetsp[ii]-&_jets[0], ymin);
2603  jetsp[ii] = jetsp[n-1];
2604  indices[ii] = indices[n-1];
2605  }
2606  }
2607 }
2609 #include<iostream>
2610 FJCORE_BEGIN_NAMESPACE // defined in fastjet/internal/base.hh
2611 using namespace std;
2612 template<> inline void ClusterSequence::_bj_set_jetinfo(
2613  EEBriefJet * const jetA, const int _jets_index) const {
2614  double E = _jets[_jets_index].E();
2615  double scale = E*E; // the default energy scale for the kt alg
2616  double p = jet_def().extra_param(); // in case we're ee_genkt
2617  switch (_jet_algorithm) {
2618  case ee_kt_algorithm:
2619  assert(_Rparam > 2.0); // force this to be true! [not best place, but works]
2620  break;
2621  case ee_genkt_algorithm:
2622  if (p <= 0 && scale < 1e-300) scale = 1e-300; // same dodgy safety as genkt
2623  scale = pow(scale,p);
2624  break;
2625  default:
2626  throw Error("Unrecognised jet algorithm");
2627  }
2628  jetA->kt2 = scale; // "kt2" might one day be renamed as "scale" or some such
2629  double norm = _jets[_jets_index].modp2();
2630  if (norm > 0) {
2631  norm = 1.0/sqrt(norm);
2632  jetA->nx = norm * _jets[_jets_index].px();
2633  jetA->ny = norm * _jets[_jets_index].py();
2634  jetA->nz = norm * _jets[_jets_index].pz();
2635  } else {
2636  jetA->nx = 0.0;
2637  jetA->ny = 0.0;
2638  jetA->nz = 1.0;
2639  }
2640  jetA->_jets_index = _jets_index;
2641  jetA->NN_dist = _R2;
2642  jetA->NN = NULL;
2643 }
2644 template<> double ClusterSequence::_bj_dist(
2645  const EEBriefJet * const jeta,
2646  const EEBriefJet * const jetb) const {
2647  double dist = 1.0
2648  - jeta->nx*jetb->nx
2649  - jeta->ny*jetb->ny
2650  - jeta->nz*jetb->nz;
2651  dist *= 2; // distance is _2_*min(Ei^2,Ej^2)*(1-cos theta)
2652  return dist;
2653 }
2655  _simple_N2_cluster<BriefJet>();
2656 }
2658  _simple_N2_cluster<EEBriefJet>();
2659 }
2661 #include <iostream>
2662 FJCORE_BEGIN_NAMESPACE // defined in fastjet/internal/base.hh
2663 using namespace std;
2665  if (_associated_cs != NULL
2666  && _associated_cs->will_delete_self_when_unused()) {
2667  _associated_cs->signal_imminent_self_deletion();
2668  delete _associated_cs;
2669  }
2670 }
2672  return (_associated_cs != NULL);
2673 }
2675  return _associated_cs;
2676 }
2677 const ClusterSequence * ClusterSequenceStructure::validated_cs() const {
2678  if (!_associated_cs)
2679  throw Error("you requested information about the internal structure of a jet, but its associated ClusterSequence has gone out of scope.");
2680  return _associated_cs;
2681 }
2682 bool ClusterSequenceStructure::has_partner(const PseudoJet &reference, PseudoJet &partner) const{
2683  return validated_cs()->has_partner(reference, partner);
2684 }
2685 bool ClusterSequenceStructure::has_child(const PseudoJet &reference, PseudoJet &child) const{
2686  return validated_cs()->has_child(reference, child);
2687 }
2688 bool ClusterSequenceStructure::has_parents(const PseudoJet &reference, PseudoJet &parent1, PseudoJet &parent2) const{
2689  return validated_cs()->has_parents(reference, parent1, parent2);
2690 }
2691 bool ClusterSequenceStructure::object_in_jet(const PseudoJet &reference, const PseudoJet &jet) const{
2692  if ((!has_associated_cluster_sequence()) || (!jet.has_associated_cluster_sequence()))
2693  throw Error("you requested information about the internal structure of a jet, but it is not associated with a ClusterSequence or its associated ClusterSequence has gone out of scope.");
2694  if (reference.associated_cluster_sequence() != jet.associated_cluster_sequence()) return false;
2695  return validated_cs()->object_in_jet(reference, jet);
2696 }
2698  if (!has_associated_cluster_sequence())
2699  throw Error("you requested information about the internal structure of a jet, but it is not associated with a ClusterSequence or its associated ClusterSequence has gone out of scope.");
2700  return true;
2701 }
2702 vector<PseudoJet> ClusterSequenceStructure::constituents(const PseudoJet &reference) const{
2703  return validated_cs()->constituents(reference);
2704 }
2706  if (!has_associated_cluster_sequence())
2707  throw Error("you requested information about the internal structure of a jet, but it is not associated with a ClusterSequence or its associated ClusterSequence has gone out of scope.");
2708  return true;
2709 }
2710 std::vector<PseudoJet> ClusterSequenceStructure::exclusive_subjets (const PseudoJet &reference, const double & dcut) const {
2711  return validated_cs()->exclusive_subjets(reference, dcut);
2712 }
2713 int ClusterSequenceStructure::n_exclusive_subjets(const PseudoJet &reference, const double & dcut) const {
2714  return validated_cs()->n_exclusive_subjets(reference, dcut);
2715 }
2716 std::vector<PseudoJet> ClusterSequenceStructure::exclusive_subjets_up_to (const PseudoJet &reference, int nsub) const {
2717  return validated_cs()->exclusive_subjets_up_to(reference, nsub);
2718 }
2719 double ClusterSequenceStructure::exclusive_subdmerge(const PseudoJet &reference, int nsub) const {
2720  return validated_cs()->exclusive_subdmerge(reference, nsub);
2721 }
2722 double ClusterSequenceStructure::exclusive_subdmerge_max(const PseudoJet &reference, int nsub) const {
2723  return validated_cs()->exclusive_subdmerge_max(reference, nsub);
2724 }
2725 bool ClusterSequenceStructure::has_pieces(const PseudoJet &reference) const{
2726  PseudoJet dummy1, dummy2;
2727  return has_parents(reference, dummy1, dummy2);
2728 }
2729 vector<PseudoJet> ClusterSequenceStructure::pieces(const PseudoJet &reference) const{
2730  PseudoJet j1, j2;
2731  vector<PseudoJet> res;
2732  if (has_parents(reference, j1, j2)){
2733  res.push_back(j1);
2734  res.push_back(j2);
2735  }
2736  return res;
2737 }
2739 #include<iostream>
2740 #include<vector>
2741 #include<cmath>
2742 #include<algorithm>
2743 FJCORE_BEGIN_NAMESPACE // defined in fastjet/internal/base.hh
2744 using namespace std;
2746  Tile * tile = & _tiles[jet->tile_index];
2747  if (jet->previous == NULL) {
2748  tile->head = jet->next;
2749  } else {
2750  jet->previous->next = jet->next;
2751  }
2752  if (jet->next != NULL) {
2753  jet->next->previous = jet->previous;
2754  }
2755 }
2757  double default_size = max(0.1,_Rparam);
2758  _tile_size_eta = default_size;
2759  _n_tiles_phi = max(3,int(floor(twopi/default_size)));
2760  _tile_size_phi = twopi / _n_tiles_phi; // >= _Rparam and fits in 2pi
2761  TilingExtent tiling_analysis(*this);
2762  _tiles_eta_min = tiling_analysis.minrap();
2763  _tiles_eta_max = tiling_analysis.maxrap();
2769  for (int ieta = _tiles_ieta_min; ieta <= _tiles_ieta_max; ieta++) {
2770  for (int iphi = 0; iphi < _n_tiles_phi; iphi++) {
2771  Tile * tile = & _tiles[_tile_index(ieta,iphi)];
2772  tile->head = NULL; // first element of tiles points to itself
2773  tile->begin_tiles[0] = tile;
2774  Tile ** pptile = & (tile->begin_tiles[0]);
2775  pptile++;
2776  tile->surrounding_tiles = pptile;
2777  if (ieta > _tiles_ieta_min) {
2778  // with the itile subroutine, we can safely run tiles from
2779  // idphi=-1 to idphi=+1, because it takes care of
2780  // negative and positive boundaries
2781  for (int idphi = -1; idphi <=+1; idphi++) {
2782  *pptile = & _tiles[_tile_index(ieta-1,iphi+idphi)];
2783  pptile++;
2784  }
2785  }
2786  *pptile = & _tiles[_tile_index(ieta,iphi-1)];
2787  pptile++;
2788  tile->RH_tiles = pptile;
2789  *pptile = & _tiles[_tile_index(ieta,iphi+1)];
2790  pptile++;
2791  if (ieta < _tiles_ieta_max) {
2792  for (int idphi = -1; idphi <= +1; idphi++) {
2793  *pptile = & _tiles[_tile_index(ieta+1,iphi+idphi)];
2794  pptile++;
2795  }
2796  }
2797  tile->end_tiles = pptile;
2798  tile->tagged = false;
2799  }
2800  }
2801 }
2802 int ClusterSequence::_tile_index(const double eta, const double phi) const {
2803  int ieta, iphi;
2804  if (eta <= _tiles_eta_min) {ieta = 0;}
2805  else if (eta >= _tiles_eta_max) {ieta = _tiles_ieta_max-_tiles_ieta_min;}
2806  else {
2807  ieta = int(((eta - _tiles_eta_min) / _tile_size_eta));
2808  if (ieta > _tiles_ieta_max-_tiles_ieta_min) {
2810  }
2811  iphi = int((phi+twopi)/_tile_size_phi) % _n_tiles_phi;
2812  return (iphi + ieta * _n_tiles_phi);
2813 }
2815  const int _jets_index) {
2816  _bj_set_jetinfo<>(jet, _jets_index);
2817  jet->tile_index = _tile_index(jet->eta, jet->phi);
2818  Tile * tile = &_tiles[jet->tile_index];
2819  jet->previous = NULL;
2820  jet->next = tile->head;
2821  if (jet->next != NULL) {jet->next->previous = jet;}
2822  tile->head = jet;
2823 }
2824 void ClusterSequence::_print_tiles(TiledJet * briefjets ) const {
2825  for (vector<Tile>::const_iterator tile = _tiles.begin();
2826  tile < _tiles.end(); tile++) {
2827  cout << "Tile " << tile - _tiles.begin()<<" = ";
2828  vector<int> list;
2829  for (TiledJet * jetI = tile->head; jetI != NULL; jetI = jetI->next) {
2830  list.push_back(jetI-briefjets);
2831  }
2832  sort(list.begin(),list.end());
2833  for (unsigned int i = 0; i < list.size(); i++) {cout <<" "<<list[i];}
2834  cout <<"\n";
2835  }
2836 }
2838  vector<int> & tile_union, int & n_near_tiles) const {
2839  for (Tile * const * near_tile = _tiles[tile_index].begin_tiles;
2840  near_tile != _tiles[tile_index].end_tiles; near_tile++){
2841  tile_union[n_near_tiles] = *near_tile - & _tiles[0];
2842  n_near_tiles++;
2843  }
2844 }
2846  const int tile_index,
2847  vector<int> & tile_union, int & n_near_tiles) {
2848  for (Tile ** near_tile = _tiles[tile_index].begin_tiles;
2849  near_tile != _tiles[tile_index].end_tiles; near_tile++){
2850  if (! (*near_tile)->tagged) {
2851  (*near_tile)->tagged = true;
2852  tile_union[n_near_tiles] = *near_tile - & _tiles[0];
2853  n_near_tiles++;
2854  }
2855  }
2856 }
2859  int n = _jets.size();
2860  TiledJet * briefjets = new TiledJet[n];
2861  TiledJet * jetA = briefjets, * jetB;
2862  TiledJet oldB;
2863  oldB.tile_index=0; // prevents a gcc warning
2864  vector<int> tile_union(3*n_tile_neighbours);
2865  for (int i = 0; i< n; i++) {
2866  _tj_set_jetinfo(jetA, i);
2867  jetA++; // move on to next entry of briefjets
2868  }
2869  TiledJet * tail = jetA; // a semaphore for the end of briefjets
2870  TiledJet * head = briefjets; // a nicer way of naming start
2871  vector<Tile>::const_iterator tile;
2872  for (tile = _tiles.begin(); tile != _tiles.end(); tile++) {
2873  for (jetA = tile->head; jetA != NULL; jetA = jetA->next) {
2874  for (jetB = tile->head; jetB != jetA; jetB = jetB->next) {
2875  double dist = _bj_dist(jetA,jetB);
2876  if (dist < jetA->NN_dist) {jetA->NN_dist = dist; jetA->NN = jetB;}
2877  if (dist < jetB->NN_dist) {jetB->NN_dist = dist; jetB->NN = jetA;}
2878  }
2879  }
2880  for (Tile ** RTile = tile->RH_tiles; RTile != tile->end_tiles; RTile++) {
2881  for (jetA = tile->head; jetA != NULL; jetA = jetA->next) {
2882  for (jetB = (*RTile)->head; jetB != NULL; jetB = jetB->next) {
2883  double dist = _bj_dist(jetA,jetB);
2884  if (dist < jetA->NN_dist) {jetA->NN_dist = dist; jetA->NN = jetB;}
2885  if (dist < jetB->NN_dist) {jetB->NN_dist = dist; jetB->NN = jetA;}
2886  }
2887  }
2888  }
2889  }
2890  double * diJ = new double[n];
2891  jetA = head;
2892  for (int i = 0; i < n; i++) {
2893  diJ[i] = _bj_diJ(jetA);
2894  jetA++; // have jetA follow i
2895  }
2896  int history_location = n-1;
2897  while (tail != head) {
2898  double diJ_min = diJ[0];
2899  int diJ_min_jet = 0;
2900  for (int i = 1; i < n; i++) {
2901  if (diJ[i] < diJ_min) {diJ_min_jet = i; diJ_min = diJ[i];}
2902  }
2903  history_location++;
2904  jetA = & briefjets[diJ_min_jet];
2905  jetB = jetA->NN;
2906  diJ_min *= _invR2;
2907  if (jetB != NULL) {
2908  if (jetA < jetB) {std::swap(jetA,jetB);}
2909  int nn; // new jet index
2910  _do_ij_recombination_step(jetA->_jets_index, jetB->_jets_index, diJ_min, nn);
2911  _bj_remove_from_tiles(jetA);
2912  oldB = * jetB; // take a copy because we will need it...
2913  _bj_remove_from_tiles(jetB);
2914  _tj_set_jetinfo(jetB, nn); // also registers the jet in the tiling
2915  } else {
2916  _do_iB_recombination_step(jetA->_jets_index, diJ_min);
2917  _bj_remove_from_tiles(jetA);
2918  }
2919  int n_near_tiles = 0;
2920  _add_neighbours_to_tile_union(jetA->tile_index, tile_union, n_near_tiles);
2921  if (jetB != NULL) {
2922  bool sort_it = false;
2923  if (jetB->tile_index != jetA->tile_index) {
2924  sort_it = true;
2925  _add_neighbours_to_tile_union(jetB->tile_index,tile_union,n_near_tiles);
2926  }
2927  if (oldB.tile_index != jetA->tile_index &&
2928  oldB.tile_index != jetB->tile_index) {
2929  sort_it = true;
2930  _add_neighbours_to_tile_union(oldB.tile_index,tile_union,n_near_tiles);
2931  }
2932  if (sort_it) {
2933  // sort the tiles before then compressing the list
2934  sort(tile_union.begin(), tile_union.begin()+n_near_tiles);
2935  // and now condense the list
2936  int nnn = 1;
2937  for (int i = 1; i < n_near_tiles; i++) {
2938  if (tile_union[i] != tile_union[nnn-1]) {
2939  tile_union[nnn] = tile_union[i];
2940  nnn++;
2941  }
2942  }
2943  n_near_tiles = nnn;
2944  }
2945  }
2946  tail--; n--;
2947  if (jetA == tail) {
2948  } else {
2949  *jetA = *tail;
2950  diJ[jetA - head] = diJ[tail-head];
2951  if (jetA->previous == NULL) {
2952  _tiles[jetA->tile_index].head = jetA;
2953  } else {
2954  jetA->previous->next = jetA;
2955  }
2956  if (jetA->next != NULL) {jetA->next->previous = jetA;}
2957  }
2958  for (int itile = 0; itile < n_near_tiles; itile++) {
2959  Tile * tile_ptr = &_tiles[tile_union[itile]];
2960  for (TiledJet * jetI = tile_ptr->head; jetI != NULL; jetI = jetI->next) {
2961  // see if jetI had jetA or jetB as a NN -- if so recalculate the NN
2962  if (jetI->NN == jetA || (jetI->NN == jetB && jetB != NULL)) {
2963  jetI->NN_dist = _R2;
2964  jetI->NN = NULL;
2965  // now go over tiles that are neighbours of I (include own tile)
2966  for (Tile ** near_tile = tile_ptr->begin_tiles;
2967  near_tile != tile_ptr->end_tiles; near_tile++) {
2968  // and then over the contents of that tile
2969  for (TiledJet * jetJ = (*near_tile)->head;
2970  jetJ != NULL; jetJ = jetJ->next) {
2971  double dist = _bj_dist(jetI,jetJ);
2972  if (dist < jetI->NN_dist && jetJ != jetI) {
2973  jetI->NN_dist = dist; jetI->NN = jetJ;
2974  }
2975  }
2976  }
2977  diJ[jetI-head] = _bj_diJ(jetI); // update diJ
2978  }
2979  // check whether new jetB is closer than jetI's current NN and
2980  // if need to update things
2981  if (jetB != NULL) {
2982  double dist = _bj_dist(jetI,jetB);
2983  if (dist < jetI->NN_dist) {
2984  if (jetI != jetB) {
2985  jetI->NN_dist = dist;
2986  jetI->NN = jetB;
2987  diJ[jetI-head] = _bj_diJ(jetI); // update diJ...
2988  }
2989  }
2990  if (dist < jetB->NN_dist) {
2991  if (jetI != jetB) {
2992  jetB->NN_dist = dist;
2993  jetB->NN = jetI;}
2994  }
2995  }
2996  }
2997  }
2998  if (jetB != NULL) {diJ[jetB-head] = _bj_diJ(jetB);}
2999  for (Tile ** near_tile = _tiles[tail->tile_index].begin_tiles;
3000  near_tile!= _tiles[tail->tile_index].end_tiles; near_tile++){
3001  for (TiledJet * jetJ = (*near_tile)->head;
3002  jetJ != NULL; jetJ = jetJ->next) {
3003  if (jetJ->NN == tail) {jetJ->NN = jetA;}
3004  }
3005  }
3006  if (jetB != NULL) {diJ[jetB-head] = _bj_diJ(jetB);}
3007  }
3008  delete[] diJ;
3009  delete[] briefjets;
3010 }
3013  int n = _jets.size();
3014  TiledJet * briefjets = new TiledJet[n];
3015  TiledJet * jetA = briefjets, * jetB;
3016  TiledJet oldB;
3017  oldB.tile_index=0; // prevents a gcc warning
3018  vector<int> tile_union(3*n_tile_neighbours);
3019  for (int i = 0; i< n; i++) {
3020  _tj_set_jetinfo(jetA, i);
3021  jetA++; // move on to next entry of briefjets
3022  }
3023  TiledJet * head = briefjets; // a nicer way of naming start
3024  vector<Tile>::const_iterator tile;
3025  for (tile = _tiles.begin(); tile != _tiles.end(); tile++) {
3026  for (jetA = tile->head; jetA != NULL; jetA = jetA->next) {
3027  for (jetB = tile->head; jetB != jetA; jetB = jetB->next) {
3028  double dist = _bj_dist(jetA,jetB);
3029  if (dist < jetA->NN_dist) {jetA->NN_dist = dist; jetA->NN = jetB;}
3030  if (dist < jetB->NN_dist) {jetB->NN_dist = dist; jetB->NN = jetA;}
3031  }
3032  }
3033  for (Tile ** RTile = tile->RH_tiles; RTile != tile->end_tiles; RTile++) {
3034  for (jetA = tile->head; jetA != NULL; jetA = jetA->next) {
3035  for (jetB = (*RTile)->head; jetB != NULL; jetB = jetB->next) {
3036  double dist = _bj_dist(jetA,jetB);
3037  if (dist < jetA->NN_dist) {jetA->NN_dist = dist; jetA->NN = jetB;}
3038  if (dist < jetB->NN_dist) {jetB->NN_dist = dist; jetB->NN = jetA;}
3039  }
3040  }
3041  }
3042  }
3043  struct diJ_plus_link {
3044  double diJ; // the distance
3045  TiledJet * jet; // the jet (i) for which we've found this distance
3046  };
3047  diJ_plus_link * diJ = new diJ_plus_link[n];
3048  jetA = head;
3049  for (int i = 0; i < n; i++) {
3050  diJ[i].diJ = _bj_diJ(jetA); // kt distance * R^2
3051  diJ[i].jet = jetA; // our compact diJ table will not be in
3052  jetA->diJ_posn = i; // one-to-one corresp. with non-compact jets,
3053  jetA++; // have jetA follow i
3054  }
3055  int history_location = n-1;
3056  while (n > 0) {
3057  diJ_plus_link * best, *stop; // pointers a bit faster than indices
3058  double diJ_min = diJ[0].diJ; // initialise the best one here.
3059  best = diJ; // and here
3060  stop = diJ+n;
3061  for (diJ_plus_link * here = diJ+1; here != stop; here++) {
3062  if (here->diJ < diJ_min) {best = here; diJ_min = here->diJ;}
3063  }
3064  history_location++;
3065  jetA = best->jet;
3066  jetB = jetA->NN;
3067  diJ_min *= _invR2;
3068  if (jetB != NULL) {
3069  if (jetA < jetB) {std::swap(jetA,jetB);}
3070  int nn; // new jet index
3071  _do_ij_recombination_step(jetA->_jets_index, jetB->_jets_index, diJ_min, nn);
3072  _bj_remove_from_tiles(jetA);
3073  oldB = * jetB; // take a copy because we will need it...
3074  _bj_remove_from_tiles(jetB);
3075  _tj_set_jetinfo(jetB, nn); // cause jetB to become _jets[nn]
3076  } else {
3077  _do_iB_recombination_step(jetA->_jets_index, diJ_min);
3078  _bj_remove_from_tiles(jetA);
3079  }
3080  int n_near_tiles = 0;
3082  tile_union, n_near_tiles);
3083  if (jetB != NULL) {
3084  if (jetB->tile_index != jetA->tile_index) {
3085  _add_untagged_neighbours_to_tile_union(jetB->tile_index,
3086  tile_union,n_near_tiles);
3087  }
3088  if (oldB.tile_index != jetA->tile_index &&
3089  oldB.tile_index != jetB->tile_index) {
3091  tile_union,n_near_tiles);
3092  }
3093  }
3094  n--;
3095  diJ[n].jet->diJ_posn = jetA->diJ_posn;
3096  diJ[jetA->diJ_posn] = diJ[n];
3097  for (int itile = 0; itile < n_near_tiles; itile++) {
3098  Tile * tile_ptr = &_tiles[tile_union[itile]];
3099  tile_ptr->tagged = false; // reset tag, since we're done with unions
3100  for (TiledJet * jetI = tile_ptr->head; jetI != NULL; jetI = jetI->next) {
3101  // see if jetI had jetA or jetB as a NN -- if so recalculate the NN
3102  if (jetI->NN == jetA || (jetI->NN == jetB && jetB != NULL)) {
3103  jetI->NN_dist = _R2;
3104  jetI->NN = NULL;
3105  // now go over tiles that are neighbours of I (include own tile)
3106  for (Tile ** near_tile = tile_ptr->begin_tiles;
3107  near_tile != tile_ptr->end_tiles; near_tile++) {
3108  // and then over the contents of that tile
3109  for (TiledJet * jetJ = (*near_tile)->head;
3110  jetJ != NULL; jetJ = jetJ->next) {
3111  double dist = _bj_dist(jetI,jetJ);
3112  if (dist < jetI->NN_dist && jetJ != jetI) {
3113  jetI->NN_dist = dist; jetI->NN = jetJ;
3114  }
3115  }
3116  }
3117  diJ[jetI->diJ_posn].diJ = _bj_diJ(jetI); // update diJ kt-dist
3118  }
3119  // check whether new jetB is closer than jetI's current NN and
3120  // if jetI is closer than jetB's current (evolving) nearest
3121  // neighbour. Where relevant update things
3122  if (jetB != NULL) {
3123  double dist = _bj_dist(jetI,jetB);
3124  if (dist < jetI->NN_dist) {
3125  if (jetI != jetB) {
3126  jetI->NN_dist = dist;
3127  jetI->NN = jetB;
3128  diJ[jetI->diJ_posn].diJ = _bj_diJ(jetI); // update diJ...
3129  }
3130  }
3131  if (dist < jetB->NN_dist) {
3132  if (jetI != jetB) {
3133  jetB->NN_dist = dist;
3134  jetB->NN = jetI;}
3135  }
3136  }
3137  }
3138  }
3139  if (jetB != NULL) {diJ[jetB->diJ_posn].diJ = _bj_diJ(jetB);}
3140  }
3141  delete[] diJ;
3142  delete[] briefjets;
3143 }
3146  int n = _jets.size();
3147  TiledJet * briefjets = new TiledJet[n];
3148  TiledJet * jetA = briefjets, * jetB;
3149  TiledJet oldB;
3150  oldB.tile_index=0; // prevents a gcc warning
3151  vector<int> tile_union(3*n_tile_neighbours);
3152  for (int i = 0; i< n; i++) {
3153  _tj_set_jetinfo(jetA, i);
3154  jetA++; // move on to next entry of briefjets
3155  }
3156  TiledJet * head = briefjets; // a nicer way of naming start
3157  vector<Tile>::const_iterator tile;
3158  for (tile = _tiles.begin(); tile != _tiles.end(); tile++) {
3159  for (jetA = tile->head; jetA != NULL; jetA = jetA->next) {
3160  for (jetB = tile->head; jetB != jetA; jetB = jetB->next) {
3161  double dist = _bj_dist(jetA,jetB);
3162  if (dist < jetA->NN_dist) {jetA->NN_dist = dist; jetA->NN = jetB;}
3163  if (dist < jetB->NN_dist) {jetB->NN_dist = dist; jetB->NN = jetA;}
3164  }
3165  }
3166  for (Tile ** RTile = tile->RH_tiles; RTile != tile->end_tiles; RTile++) {
3167  for (jetA = tile->head; jetA != NULL; jetA = jetA->next) {
3168  for (jetB = (*RTile)->head; jetB != NULL; jetB = jetB->next) {
3169  double dist = _bj_dist(jetA,jetB);
3170  if (dist < jetA->NN_dist) {jetA->NN_dist = dist; jetA->NN = jetB;}
3171  if (dist < jetB->NN_dist) {jetB->NN_dist = dist; jetB->NN = jetA;}
3172  }
3173  }
3174  }
3175  }
3176  vector<double> diJs(n);
3177  for (int i = 0; i < n; i++) {
3178  diJs[i] = _bj_diJ(&briefjets[i]);
3179  briefjets[i].label_minheap_update_done();
3180  }
3181  MinHeap minheap(diJs);
3182  vector<TiledJet *> jets_for_minheap;
3183  jets_for_minheap.reserve(n);
3184  int history_location = n-1;
3185  while (n > 0) {
3186  double diJ_min = minheap.minval() *_invR2;
3187  jetA = head + minheap.minloc();
3188  history_location++;
3189  jetB = jetA->NN;
3190  if (jetB != NULL) {
3191  if (jetA < jetB) {std::swap(jetA,jetB);}
3192  int nn; // new jet index
3193  _do_ij_recombination_step(jetA->_jets_index, jetB->_jets_index, diJ_min, nn);
3194  _bj_remove_from_tiles(jetA);
3195  oldB = * jetB; // take a copy because we will need it...
3196  _bj_remove_from_tiles(jetB);
3197  _tj_set_jetinfo(jetB, nn); // cause jetB to become _jets[nn]
3198  } else {
3199  _do_iB_recombination_step(jetA->_jets_index, diJ_min);
3200  _bj_remove_from_tiles(jetA);
3201  }
3202  minheap.remove(jetA-head);
3203  int n_near_tiles = 0;
3205  tile_union, n_near_tiles);
3206  if (jetB != NULL) {
3207  if (jetB->tile_index != jetA->tile_index) {
3208  _add_untagged_neighbours_to_tile_union(jetB->tile_index,
3209  tile_union,n_near_tiles);
3210  }
3211  if (oldB.tile_index != jetA->tile_index &&
3212  oldB.tile_index != jetB->tile_index) {
3213  // GS: the line below generates a warning that oldB.tile_index
3214  // may be used uninitialised. However, to reach this point, we
3215  // ned jetB != NULL (see test a few lines above) and is jetB
3216  // !=NULL, one would have gone through "oldB = *jetB before
3217  // (see piece of code ~20 line above), so the index is
3218  // initialised. We do not do anything to avoid the warning to
3219  // avoid any potential speed impact.
3221  tile_union,n_near_tiles);
3222  }
3223  jetB->label_minheap_update_needed();
3224  jets_for_minheap.push_back(jetB);
3225  }
3226  for (int itile = 0; itile < n_near_tiles; itile++) {
3227  Tile * tile_ptr = &_tiles[tile_union[itile]];
3228  tile_ptr->tagged = false; // reset tag, since we're done with unions
3229  for (TiledJet * jetI = tile_ptr->head; jetI != NULL; jetI = jetI->next) {
3230  // see if jetI had jetA or jetB as a NN -- if so recalculate the NN
3231  if (jetI->NN == jetA || (jetI->NN == jetB && jetB != NULL)) {
3232  jetI->NN_dist = _R2;
3233  jetI->NN = NULL;
3234  // label jetI as needing heap action...
3235  if (!jetI->minheap_update_needed()) {
3236  jetI->label_minheap_update_needed();
3237  jets_for_minheap.push_back(jetI);}
3238  // now go over tiles that are neighbours of I (include own tile)
3239  for (Tile ** near_tile = tile_ptr->begin_tiles;
3240  near_tile != tile_ptr->end_tiles; near_tile++) {
3241  // and then over the contents of that tile
3242  for (TiledJet * jetJ = (*near_tile)->head;
3243  jetJ != NULL; jetJ = jetJ->next) {
3244  double dist = _bj_dist(jetI,jetJ);
3245  if (dist < jetI->NN_dist && jetJ != jetI) {
3246  jetI->NN_dist = dist; jetI->NN = jetJ;
3247  }
3248  }
3249  }
3250  }
3251  // check whether new jetB is closer than jetI's current NN and
3252  // if jetI is closer than jetB's current (evolving) nearest
3253  // neighbour. Where relevant update things
3254  if (jetB != NULL) {
3255  double dist = _bj_dist(jetI,jetB);
3256  if (dist < jetI->NN_dist) {
3257  if (jetI != jetB) {
3258  jetI->NN_dist = dist;
3259  jetI->NN = jetB;
3260  // label jetI as needing heap action...
3261  if (!jetI->minheap_update_needed()) {
3262  jetI->label_minheap_update_needed();
3263  jets_for_minheap.push_back(jetI);}
3264  }
3265  }
3266  if (dist < jetB->NN_dist) {
3267  if (jetI != jetB) {
3268  jetB->NN_dist = dist;
3269  jetB->NN = jetI;}
3270  }
3271  }
3272  }
3273  }
3274  while (jets_for_minheap.size() > 0) {
3275  TiledJet * jetI = jets_for_minheap.back();
3276  jets_for_minheap.pop_back();
3277  minheap.update(jetI-head, _bj_diJ(jetI));
3278  jetI->label_minheap_update_done();
3279  }
3280  n--;
3281  }
3282  delete[] briefjets;
3283 }
3285 FJCORE_BEGIN_NAMESPACE // defined in fastjet/internal/base.hh
3286 using namespace std;
3287 CompositeJetStructure::CompositeJetStructure(const std::vector<PseudoJet> & initial_pieces,
3288  const JetDefinition::Recombiner * recombiner)
3289  : _pieces(initial_pieces){
3290  if (recombiner){}; // ugly trick to prevent a gcc warning
3291  _area_4vector_ptr = 0;
3292 }
3294  string str = "Composite PseudoJet";
3295  return str;
3296 }
3298  return _pieces.size()!=0;
3299 }
3300 std::vector<PseudoJet> CompositeJetStructure::constituents(const PseudoJet & /*jet*/) const{
3301  vector<PseudoJet> all_constituents;
3302  for (unsigned i = 0; i < _pieces.size(); i++) {
3303  if (_pieces[i].has_constituents()){
3304  vector<PseudoJet> constits = _pieces[i].constituents();
3305  copy(constits.begin(), constits.end(), back_inserter(all_constituents));
3306  } else {
3307  all_constituents.push_back(_pieces[i]);
3308  }
3309  }
3310  return all_constituents;
3311 }
3312 std::vector<PseudoJet> CompositeJetStructure::pieces(const PseudoJet & /*jet*/) const{
3313  return _pieces;
3314 }
3315 FJCORE_END_NAMESPACE // defined in fastjet/internal/base.hh
3316 #include <sstream>
3317 FJCORE_BEGIN_NAMESPACE // defined in fastjet/internal/base.hh
3318 using namespace std;
3319 bool Error::_print_errors = true;
3320 bool Error::_print_backtrace = false;
3321 ostream * Error::_default_ostr = & cerr;
3322 #if (!defined(FJCORE_HAVE_EXECINFO_H)) || defined(__FJCORE__)
3324 #endif
3325 Error::Error(const std::string & message_in) {
3326  _message = message_in;
3327  if (_print_errors && _default_ostr){
3328  ostringstream oss;
3329  oss << "fjcore::Error: "<< message_in << endl;
3330  *_default_ostr << oss.str();
3331  _default_ostr->flush();
3332  }
3333 }
3334 void Error::set_print_backtrace(bool enabled) {
3335 #if (!defined(FJCORE_HAVE_EXECINFO_H)) || defined(__FJCORE__)
3336  if (enabled) {
3337  _execinfo_undefined.warn("Error::set_print_backtrace(true) will not work with this build of FastJet");
3338  }
3339 #endif
3340  _print_backtrace = enabled;
3341 }
3343 #include <string>
3344 #include <sstream>
3345 using namespace std;
3348 #include<sstream>
3349 FJCORE_BEGIN_NAMESPACE // defined in fastjet/internal/base.hh
3350 using namespace std;
3351 const double JetDefinition::max_allowable_R = 1000.0;
3352 JetDefinition::JetDefinition(JetAlgorithm jet_algorithm_in,
3353  double R_in,
3354  RecombinationScheme recomb_scheme_in,
3355  Strategy strategy_in,
3356  int nparameters) :
3357  _jet_algorithm(jet_algorithm_in), _Rparam(R_in), _strategy(strategy_in) {
3358  if (_jet_algorithm == ee_kt_algorithm) {
3359  _Rparam = 4.0; // introduce a fictional R that ensures that
3360  } else {
3361  if (R_in > max_allowable_R) {
3362  ostringstream oss;
3363  oss << "Requested R = " << R_in << " for jet definition is larger than max_allowable_R = " << max_allowable_R;
3364  throw Error(oss.str());
3365  }
3366  }
3367  unsigned int nparameters_expected = n_parameters_for_algorithm(jet_algorithm_in);
3368  if (nparameters != (int) nparameters_expected){
3369  ostringstream oss;
3370  oss << "The jet algorithm you requested ("
3371  << jet_algorithm_in << ") should be constructed with " << nparameters_expected
3372  << " parameter(s) but was called with " << nparameters << " parameter(s)\n";
3373  throw Error(oss.str());
3374  }
3375  assert (_strategy != plugin_strategy);
3376  _plugin = NULL;
3377  set_recombination_scheme(recomb_scheme_in);
3378  set_extra_param(0.0); // make sure it's defined
3379 }
3381  if (jet_algorithm() == plugin_algorithm) {
3382  return plugin()->is_spherical();
3383  } else {
3384  return (jet_algorithm() == ee_kt_algorithm || // as of 2013-02-14, the two
3385  jet_algorithm() == ee_genkt_algorithm // native spherical algorithms
3386  );
3387  }
3388 }
3390  ostringstream name;
3391  name << description_no_recombiner();
3392  if ((jet_algorithm() == plugin_algorithm) || (jet_algorithm() == undefined_jet_algorithm)){
3393  return name.str();
3394  }
3395  if (n_parameters_for_algorithm(jet_algorithm()) == 0)
3396  name << " with ";
3397  else
3398  name << " and ";
3399  name << recombiner()->description();
3400  return name.str();
3401 }
3403  ostringstream name;
3404  if (jet_algorithm() == plugin_algorithm) {
3405  return plugin()->description();
3406  } else if (jet_algorithm() == undefined_jet_algorithm) {
3407  return "uninitialised JetDefinition (jet_algorithm=undefined_jet_algorithm)" ;
3408  }
3409  name << algorithm_description(jet_algorithm());
3410  switch (n_parameters_for_algorithm(jet_algorithm())){
3411  case 0: name << " (NB: no R)"; break;
3412  case 1: name << " with R = " << R(); break; // the parameter is always R
3413  case 2:
3414  name << " with R = " << R();
3415  if (jet_algorithm() == cambridge_for_passive_algorithm){
3416  name << "and a special hack whereby particles with kt < "
3417  << extra_param() << "are treated as passive ghosts";
3418  } else {
3419  name << ", p = " << extra_param();
3420  }
3421  };
3422  return name.str();
3423 }
3425  ostringstream name;
3426  switch (jet_alg){
3427  case plugin_algorithm: return "plugin algorithm";
3428  case kt_algorithm: return "Longitudinally invariant kt algorithm";
3429  case cambridge_algorithm: return "Longitudinally invariant Cambridge/Aachen algorithm";
3430  case antikt_algorithm: return "Longitudinally invariant anti-kt algorithm";
3431  case genkt_algorithm: return "Longitudinally invariant generalised kt algorithm";
3432  case cambridge_for_passive_algorithm: return "Longitudinally invariant Cambridge/Aachen algorithm";
3433  case ee_kt_algorithm: return "e+e- kt (Durham) algorithm (NB: no R)";
3434  case ee_genkt_algorithm: return "e+e- generalised kt algorithm";
3435  case undefined_jet_algorithm: return "undefined jet algorithm";
3436  default:
3437  throw Error("JetDefinition::algorithm_description(): unrecognized jet_algorithm");
3438  };
3439 }
3441  switch (jet_alg) {
3442  case ee_kt_algorithm: return 0;
3443  case genkt_algorithm:
3444  case ee_genkt_algorithm: return 2;
3445  default: return 1;
3446  };
3447 }
3449  RecombinationScheme recomb_scheme) {
3450  _default_recombiner = JetDefinition::DefaultRecombiner(recomb_scheme);
3451  if (_shared_recombiner) _shared_recombiner.reset();
3452  _recombiner = 0;
3453 }
3454 void JetDefinition::set_recombiner(const JetDefinition &other_jet_def){
3455  assert(other_jet_def._recombiner ||
3456  other_jet_def.recombination_scheme() != external_scheme);
3457  if (other_jet_def._recombiner == 0){
3458  set_recombination_scheme(other_jet_def.recombination_scheme());
3459  return;
3460  }
3461  _recombiner = other_jet_def._recombiner;
3462  _default_recombiner = DefaultRecombiner(external_scheme);
3463  _shared_recombiner.reset(other_jet_def._shared_recombiner);
3464 }
3465 bool JetDefinition::has_same_recombiner(const JetDefinition &other_jd) const{
3466  const RecombinationScheme & scheme = recombination_scheme();
3467  if (other_jd.recombination_scheme() != scheme) return false;
3468  return (scheme != external_scheme)
3469  || (recombiner() == other_jd.recombiner());
3470 }
3472  if (_recombiner == 0){
3473  throw Error("tried to call JetDefinition::delete_recombiner_when_unused() for a JetDefinition without a user-defined recombination scheme");
3474  } else if (_shared_recombiner.get()) {
3475  throw Error("Error in JetDefinition::delete_recombiner_when_unused: the recombiner is already scheduled for deletion when unused (or was already set as shared)");
3476  }
3477  _shared_recombiner.reset(_recombiner);
3478 }
3480  if (_plugin == 0){
3481  throw Error("tried to call JetDefinition::delete_plugin_when_unused() for a JetDefinition without a plugin");
3482  }
3483  _plugin_shared.reset(_plugin);
3484 }
3486  switch(_recomb_scheme) {
3487  case E_scheme:
3488  return "E scheme recombination";
3489  case pt_scheme:
3490  return "pt scheme recombination";
3491  case pt2_scheme:
3492  return "pt2 scheme recombination";
3493  case Et_scheme:
3494  return "Et scheme recombination";
3495  case Et2_scheme:
3496  return "Et2 scheme recombination";
3497  case BIpt_scheme:
3498  return "boost-invariant pt scheme recombination";
3499  case BIpt2_scheme:
3500  return "boost-invariant pt2 scheme recombination";
3501  case WTA_pt_scheme:
3502  return "pt-ordered Winner-Takes-All recombination";
3503  case WTA_modp_scheme:
3504  return "|3-momentum|-ordered Winner-Takes-All recombination";
3505  default:
3506  ostringstream err;
3507  err << "DefaultRecombiner: unrecognized recombination scheme "
3508  << _recomb_scheme;
3509  throw Error(err.str());
3510  }
3511 }
3513  const PseudoJet & pa, const PseudoJet & pb,
3514  PseudoJet & pab) const {
3515  double weighta, weightb;
3516  switch(_recomb_scheme) {
3517  case E_scheme:
3518  pab.reset(pa.px()+pb.px(),
3519  pa.py()+pb.py(),
3520  pa.pz()+pb.pz(),
3521  pa.E ()+pb.E ());
3522  return;
3523  case pt_scheme:
3524  case Et_scheme:
3525  case BIpt_scheme:
3526  weighta = pa.perp();
3527  weightb = pb.perp();
3528  break;
3529  case pt2_scheme:
3530  case Et2_scheme:
3531  case BIpt2_scheme:
3532  weighta = pa.perp2();
3533  weightb = pb.perp2();
3534  break;
3535  case WTA_pt_scheme:{
3536  const PseudoJet & phard = (pa.pt2() >= pb.pt2()) ? pa : pb;
3537  pab.reset_PtYPhiM(pa.pt()+pb.pt(),
3538  phard.rap(), phard.phi(), phard.m());
3539  return;}
3540  case WTA_modp_scheme:{
3541  bool a_hardest = (pa.modp2() >= pb.modp2());
3542  const PseudoJet & phard = a_hardest ? pa : pb;
3543  const PseudoJet & psoft = a_hardest ? pb : pa;
3544  double modp_hard = phard.modp();
3545  double modp_ab = modp_hard + psoft.modp();
3546  if (phard.modp2()==0.0){
3547  pab.reset(0.0, 0.0, 0.0, phard.m());
3548  } else {
3549  double scale = modp_ab/modp_hard;
3550  pab.reset(phard.px()*scale, phard.py()*scale, phard.pz()*scale,
3551  sqrt(modp_ab*modp_ab + phard.m2()));
3552  }
3553  return;}
3554  default:
3555  ostringstream err;
3556  err << "DefaultRecombiner: unrecognized recombination scheme "
3557  << _recomb_scheme;
3558  throw Error(err.str());
3559  }
3560  double perp_ab = pa.perp() + pb.perp();
3561  if (perp_ab != 0.0) { // weights also non-zero...
3562  double y_ab = (weighta * pa.rap() + weightb * pb.rap())/(weighta+weightb);
3563  double phi_a = pa.phi(), phi_b = pb.phi();
3564  if (phi_a - phi_b > pi) phi_b += twopi;
3565  if (phi_a - phi_b < -pi) phi_b -= twopi;
3566  double phi_ab = (weighta * phi_a + weightb * phi_b)/(weighta+weightb);
3567  pab.reset_PtYPhiM(perp_ab,y_ab,phi_ab);
3568  } else { // weights are zero
3569  pab.reset(0.0, 0.0, 0.0, 0.0);
3570  }
3571 }
3573  switch(_recomb_scheme) {
3574  case E_scheme:
3575  case BIpt_scheme:
3576  case BIpt2_scheme:
3577  case WTA_pt_scheme:
3578  case WTA_modp_scheme:
3579  break;
3580  case pt_scheme:
3581  case pt2_scheme:
3582  {
3583  double newE = sqrt(p.perp2()+p.pz()*p.pz());
3584  p.reset_momentum(p.px(), p.py(), p.pz(), newE);
3585  }
3586  break;
3587  case Et_scheme:
3588  case Et2_scheme:
3589  {
3590  double rescale = p.E()/sqrt(p.perp2()+p.pz()*p.pz());
3591  p.reset_momentum(rescale*p.px(), rescale*p.py(), rescale*p.pz(), p.E());
3592  }
3593  break;
3594  default:
3595  ostringstream err;
3596  err << "DefaultRecombiner: unrecognized recombination scheme "
3597  << _recomb_scheme;
3598  throw Error(err.str());
3599  }
3600 }
3602  throw Error("set_ghost_separation_scale not supported");
3603 }
3604 PseudoJet join(const vector<PseudoJet> & pieces, const JetDefinition::Recombiner & recombiner){
3605  PseudoJet result; // automatically initialised to 0
3606  if (pieces.size()>0){
3607  result = pieces[0];
3608  for (unsigned int i=1; i<pieces.size(); i++)
3609  recombiner.plus_equal(result, pieces[i]);
3610  }
3611  CompositeJetStructure *cj_struct = new CompositeJetStructure(pieces, &recombiner);
3612  result.set_structure_shared_ptr(SharedPtr<PseudoJetStructureBase>(cj_struct));
3613  return result;
3614 }
3615 PseudoJet join(const PseudoJet & j1,
3617  return join(vector<PseudoJet>(1,j1), recombiner);
3618 }
3619 PseudoJet join(const PseudoJet & j1, const PseudoJet & j2,
3621  vector<PseudoJet> pieces;
3622  pieces.push_back(j1);
3623  pieces.push_back(j2);
3624  return join(pieces, recombiner);
3625 }
3626 PseudoJet join(const PseudoJet & j1, const PseudoJet & j2, const PseudoJet & j3,
3628  vector<PseudoJet> pieces;
3629  pieces.push_back(j1);
3630  pieces.push_back(j2);
3631  pieces.push_back(j3);
3632  return join(pieces, recombiner);
3633 }
3634 PseudoJet join(const PseudoJet & j1, const PseudoJet & j2, const PseudoJet & j3, const PseudoJet & j4,
3636  vector<PseudoJet> pieces;
3637  pieces.push_back(j1);
3638  pieces.push_back(j2);
3639  pieces.push_back(j3);
3640  pieces.push_back(j4);
3641  return join(pieces, recombiner);
3642 }
3644 #include <sstream>
3645 #include <limits>
3646 using namespace std;
3648 ostream * LimitedWarning::_default_ostr = &cerr;
3649 std::list< LimitedWarning::Summary > LimitedWarning::_global_warnings_summary;
3651 void LimitedWarning::warn(const char * warning, std::ostream * ostr) {
3652  if (_this_warning_summary == 0) {
3653  _global_warnings_summary.push_back(Summary(warning, 0));
3654  _this_warning_summary = & (_global_warnings_summary.back());
3655  }
3656  if (_n_warn_so_far < _max_warn) {
3657  ostringstream warnstr;
3658  warnstr << "WARNING from FastJet: ";
3659  warnstr << warning;
3660  _n_warn_so_far++;
3661  if (_n_warn_so_far == _max_warn) warnstr << " (LAST SUCH WARNING)";
3662  warnstr << std::endl;
3663  if (ostr) {
3664  (*ostr) << warnstr.str();
3665  ostr->flush(); // get something written to file even if the program aborts
3666  }
3667  }
3668  if (_this_warning_summary->second < numeric_limits<unsigned>::max()) {
3669  _this_warning_summary->second++;
3670  }
3671 }
3673  ostringstream str;
3674  for (list<Summary>::const_iterator it = _global_warnings_summary.begin();
3675  it != _global_warnings_summary.end(); it++) {
3676  str << it->second << " times: " << it->first << endl;
3677  }
3678  return str.str();
3679 }
3681 #include<iostream>
3682 #include<cmath>
3683 #include<limits>
3684 FJCORE_BEGIN_NAMESPACE // defined in fastjet/internal/base.hh
3685 using namespace std;
3686 void MinHeap::initialise(const std::vector<double> & values){
3687  for (unsigned i = values.size(); i < _heap.size(); i++) {
3688  _heap[i].value = std::numeric_limits<double>::max();
3689  _heap[i].minloc = &(_heap[i]);
3690  }
3691  for (unsigned i = 0; i < values.size(); i++) {
3692  _heap[i].value = values[i];
3693  _heap[i].minloc = &(_heap[i]);
3694  }
3695  for (unsigned i = _heap.size()-1; i > 0; i--) {
3696  ValueLoc * parent = &(_heap[(i-1)/2]);
3697  ValueLoc * here = &(_heap[i]);
3698  if (here->minloc->value < parent->minloc->value) {
3699  parent->minloc = here->minloc;
3700  }
3701  }
3702 }
3703 void MinHeap::update(unsigned int loc, double new_value) {
3704  assert(loc < _heap.size());
3705  ValueLoc * start = &(_heap[loc]);
3706  if (start->minloc != start && !(new_value < start->minloc->value)) {
3707  start->value = new_value;
3708  return;
3709  }
3710  start->value = new_value;
3711  start->minloc = start;
3712  bool change_made = true;
3713  ValueLoc * heap_end = (&(_heap[0])) + _heap.size();
3714  while(change_made) {
3715  ValueLoc * here = &(_heap[loc]);
3716  change_made = false;
3717  if (here->minloc == start) {
3718  here->minloc = here; change_made = true;
3719  }
3720  ValueLoc * child = &(_heap[2*loc+1]);
3721  if (child < heap_end && child->minloc->value < here->minloc->value ) {
3722  here->minloc = child->minloc;
3723  change_made = true;}
3724  child++;
3725  if (child < heap_end && child->minloc->value < here->minloc->value ) {
3726  here->minloc = child->minloc;
3727  change_made = true;}
3728  if (loc == 0) {break;}
3729  loc = (loc-1)/2;
3730  }
3731 }
3733 #include<valarray>
3734 #include<iostream>
3735 #include<sstream>
3736 #include<cmath>
3737 #include<algorithm>
3738 #include <cstdarg>
3739 FJCORE_BEGIN_NAMESPACE // defined in fastjet/internal/base.hh
3740 using namespace std;
3741 PseudoJet::PseudoJet(const double px_in, const double py_in, const double pz_in, const double E_in) {
3742  _E = E_in ;
3743  _px = px_in;
3744  _py = py_in;
3745  _pz = pz_in;
3746  this->_finish_init();
3747  _reset_indices();
3748 }
3750  _kt2 = this->px()*this->px() + this->py()*this->py();
3751  _phi = pseudojet_invalid_phi;
3752  _rap = pseudojet_invalid_rap;
3753 }
3755  if (_kt2 == 0.0) {
3756  _phi = 0.0; }
3757  else {
3758  _phi = atan2(this->py(),this->px());
3759  }
3760  if (_phi < 0.0) {_phi += twopi;}
3761  if (_phi >= twopi) {_phi -= twopi;} // can happen if phi=-|eps<1e-15|?
3762  if (this->E() == abs(this->pz()) && _kt2 == 0) {
3763  double MaxRapHere = MaxRap + abs(this->pz());
3764  if (this->pz() >= 0.0) {_rap = MaxRapHere;} else {_rap = -MaxRapHere;}
3765  } else {
3766  double effective_m2 = max(0.0,m2()); // force non tachyonic mass
3767  double E_plus_pz = _E + abs(_pz); // the safer of p+, p-
3768  _rap = 0.5*log((_kt2 + effective_m2)/(E_plus_pz*E_plus_pz));
3769  if (_pz > 0) {_rap = - _rap;}
3770  }
3771 }
3772 valarray<double> PseudoJet::four_mom() const {
3773  valarray<double> mom(4);
3774  mom[0] = _px;
3775  mom[1] = _py;
3776  mom[2] = _pz;
3777  mom[3] = _E ;
3778  return mom;
3779 }
3780 double PseudoJet::operator () (int i) const {
3781  switch(i) {
3782  case X:
3783  return px();
3784  case Y:
3785  return py();
3786  case Z:
3787  return pz();
3788  case T:
3789  return e();
3790  default:
3791  ostringstream err;
3792  err << "PseudoJet subscripting: bad index (" << i << ")";
3793  throw Error(err.str());
3794  }
3795  return 0.;
3796 }
3798  if (px() == 0.0 && py() ==0.0) return MaxRap;
3799  if (pz() == 0.0) return 0.0;
3800  double theta = atan(perp()/pz());
3801  if (theta < 0) theta += pi;
3802  return -log(tan(theta/2));
3803 }
3804 PseudoJet operator+ (const PseudoJet & jet1, const PseudoJet & jet2) {
3805  return PseudoJet(jet1.px()+jet2.px(),
3806  jet1.py()+jet2.py(),
3807  jet1.pz()+jet2.pz(),
3808  jet1.E() +jet2.E() );
3809 }
3810 PseudoJet operator- (const PseudoJet & jet1, const PseudoJet & jet2) {
3811  return PseudoJet(jet1.px()-jet2.px(),
3812  jet1.py()-jet2.py(),
3813  jet1.pz()-jet2.pz(),
3814  jet1.E() -jet2.E() );
3815 }
3816 PseudoJet operator* (double coeff, const PseudoJet & jet) {
3817  jet._ensure_valid_rap_phi();
3818  PseudoJet coeff_times_jet(jet);
3819  coeff_times_jet *= coeff;
3820  return coeff_times_jet;
3821 }
3822 PseudoJet operator* (const PseudoJet & jet, double coeff) {
3823  return coeff*jet;
3824 }
3825 PseudoJet operator/ (const PseudoJet & jet, double coeff) {
3826  return (1.0/coeff)*jet;
3827 }
3828 void PseudoJet::operator*=(double coeff) {
3829  _ensure_valid_rap_phi();
3830  _px *= coeff;
3831  _py *= coeff;
3832  _pz *= coeff;
3833  _E *= coeff;
3834  _kt2*= coeff*coeff;
3835 }
3836 void PseudoJet::operator/=(double coeff) {
3837  (*this) *= 1.0/coeff;
3838 }
3839 void PseudoJet::operator+=(const PseudoJet & other_jet) {
3840  _px += other_jet._px;
3841  _py += other_jet._py;
3842  _pz += other_jet._pz;
3843  _E += other_jet._E ;
3844  _finish_init(); // we need to recalculate phi,rap,kt2
3845 }
3846 void PseudoJet::operator-=(const PseudoJet & other_jet) {
3847  _px -= other_jet._px;
3848  _py -= other_jet._py;
3849  _pz -= other_jet._pz;
3850  _E -= other_jet._E ;
3851  _finish_init(); // we need to recalculate phi,rap,kt2
3852 }
3853 bool operator==(const PseudoJet & a, const PseudoJet & b) {
3854  if (a.px() != b.px()) return false;
3855  if (a.py() != b.py()) return false;
3856  if (a.pz() != b.pz()) return false;
3857  if (a.E () != b.E ()) return false;
3858  if (a.user_index() != b.user_index()) return false;
3859  if (a.cluster_hist_index() != b.cluster_hist_index()) return false;
3860  if (a.user_info_ptr() != b.user_info_ptr()) return false;
3861  if (a.structure_ptr() != b.structure_ptr()) return false;
3862  return true;
3863 }
3864 bool operator==(const PseudoJet & jet, const double val) {
3865  if (val != 0)
3866  throw Error("comparing a PseudoJet with a non-zero constant (double) is not allowed.");
3867  return (jet.px() == 0 && jet.py() == 0 &&
3868  jet.pz() == 0 && jet.E() == 0);
3869 }
3870 PseudoJet & PseudoJet::boost(const PseudoJet & prest) {
3871  if (prest.px() == 0.0 && prest.py() == 0.0 && prest.pz() == 0.0)
3872  return *this;
3873  double m_local = prest.m();
3874  assert(m_local != 0);
3875  double pf4 = ( px()*prest.px() + py()*prest.py()
3876  + pz()*prest.pz() + E()*prest.E() )/m_local;
3877  double fn = (pf4 + E()) / (prest.E() + m_local);
3878  _px += fn*prest.px();
3879  _py += fn*prest.py();
3880  _pz += fn*prest.pz();
3881  _E = pf4;
3882  _finish_init(); // we need to recalculate phi,rap,kt2
3883  return *this;
3884 }
3885 PseudoJet & PseudoJet::unboost(const PseudoJet & prest) {
3886  if (prest.px() == 0.0 && prest.py() == 0.0 && prest.pz() == 0.0)
3887  return *this;
3888  double m_local = prest.m();
3889  assert(m_local != 0);
3890  double pf4 = ( -px()*prest.px() - py()*prest.py()
3891  - pz()*prest.pz() + E()*prest.E() )/m_local;
3892  double fn = (pf4 + E()) / (prest.E() + m_local);
3893  _px -= fn*prest.px();
3894  _py -= fn*prest.py();
3895  _pz -= fn*prest.pz();
3896  _E = pf4;
3897  _finish_init(); // we need to recalculate phi,rap,kt2
3898  return *this;
3899 }
3900 bool have_same_momentum(const PseudoJet & jeta, const PseudoJet & jetb) {
3901  return jeta.px() == jetb.px()
3902  && jeta.py() == jetb.py()
3903  && jeta.pz() == jetb.pz()
3904  && jeta.E() == jetb.E();
3905 }
3906 void PseudoJet::set_cached_rap_phi(double rap_in, double phi_in) {
3907  _rap = rap_in; _phi = phi_in;
3908  if (_phi >= twopi) _phi -= twopi;
3909  if (_phi < 0) _phi += twopi;
3910 }
3911 void PseudoJet::reset_momentum_PtYPhiM(double pt_in, double y_in, double phi_in, double m_in) {
3912  assert(phi_in < 2*twopi && phi_in > -twopi);
3913  double ptm = (m_in == 0) ? pt_in : sqrt(pt_in*pt_in+m_in*m_in);
3914  double exprap = exp(y_in);
3915  double pminus = ptm/exprap;
3916  double pplus = ptm*exprap;
3917  double px_local = pt_in*cos(phi_in);
3918  double py_local = pt_in*sin(phi_in);
3919  reset_momentum(px_local,py_local,0.5*(pplus-pminus),0.5*(pplus+pminus));
3920  set_cached_rap_phi(y_in,phi_in);
3921 }
3922 PseudoJet PtYPhiM(double pt, double y, double phi, double m) {
3923  assert(phi < 2*twopi && phi > -twopi);
3924  double ptm = (m == 0) ? pt : sqrt(pt*pt+m*m);
3925  double exprap = exp(y);
3926  double pminus = ptm/exprap;
3927  double pplus = ptm*exprap;
3928  double px = pt*cos(phi);
3929  double py = pt*sin(phi);
3930  PseudoJet mom(px,py,0.5*(pplus-pminus),0.5*(pplus+pminus));
3931  mom.set_cached_rap_phi(y,phi);
3932  return mom;
3933 }
3934 double PseudoJet::kt_distance(const PseudoJet & other) const {
3935  double distance = min(_kt2, other._kt2);
3936  double dphi = abs(phi() - other.phi());
3937  if (dphi > pi) {dphi = twopi - dphi;}
3938  double drap = rap() - other.rap();
3939  distance = distance * (dphi*dphi + drap*drap);
3940  return distance;
3941 }
3942 double PseudoJet::plain_distance(const PseudoJet & other) const {
3943  double dphi = abs(phi() - other.phi());
3944  if (dphi > pi) {dphi = twopi - dphi;}
3945  double drap = rap() - other.rap();
3946  return (dphi*dphi + drap*drap);
3947 }
3948 double PseudoJet::delta_phi_to(const PseudoJet & other) const {
3949  double dphi = other.phi() - phi();
3950  if (dphi > pi) dphi -= twopi;
3951  if (dphi < -pi) dphi += twopi;
3952  return dphi;
3953 }
3954 string PseudoJet::description() const{
3955  if (!_structure)
3956  return "standard PseudoJet (with no associated clustering information)";
3957  return _structure->description();
3958 }
3960  return (_structure) && (_structure->has_associated_cluster_sequence());
3961 }
3962 const ClusterSequence* PseudoJet::associated_cluster_sequence() const{
3963  if (! has_associated_cluster_sequence()) return NULL;
3964  return _structure->associated_cluster_sequence();
3965 }
3967  return (_structure) && (_structure->has_valid_cluster_sequence());
3968 }
3969 const ClusterSequence * PseudoJet::validated_cs() const {
3970  return validated_structure_ptr()->validated_cs();
3971 }
3973  _structure = structure_in;
3974 }
3976  return bool(_structure);
3977 }
3979  return _structure.get();
3980 }
3982  return _structure.get();
3983 }
3985  if (!_structure)
3986  throw Error("Trying to access the structure of a PseudoJet which has no associated structure");
3987  return _structure.get();
3988 }
3990  return _structure;
3991 }
3992 bool PseudoJet::has_partner(PseudoJet &partner) const{
3993  return validated_structure_ptr()->has_partner(*this, partner);
3994 }
3995 bool PseudoJet::has_child(PseudoJet &child) const{
3996  return validated_structure_ptr()->has_child(*this, child);
3997 }
3998 bool PseudoJet::has_parents(PseudoJet &parent1, PseudoJet &parent2) const{
3999  return validated_structure_ptr()->has_parents(*this, parent1, parent2);
4000 }
4001 bool PseudoJet::contains(const PseudoJet &constituent) const{
4002  return validated_structure_ptr()->object_in_jet(constituent, *this);
4003 }
4004 bool PseudoJet::is_inside(const PseudoJet &jet) const{
4005  return validated_structure_ptr()->object_in_jet(*this, jet);
4006 }
4008  return (_structure) && (_structure->has_constituents());
4009 }
4010 vector<PseudoJet> PseudoJet::constituents() const{
4011  return validated_structure_ptr()->constituents(*this);
4012 }
4014  return (_structure) && (_structure->has_exclusive_subjets());
4015 }
4016 std::vector<PseudoJet> PseudoJet::exclusive_subjets (const double dcut) const {
4017  return validated_structure_ptr()->exclusive_subjets(*this, dcut);
4018 }
4019 int PseudoJet::n_exclusive_subjets(const double dcut) const {
4020  return validated_structure_ptr()->n_exclusive_subjets(*this, dcut);
4021 }
4022 std::vector<PseudoJet> PseudoJet::exclusive_subjets_up_to (int nsub) const {
4023  return validated_structure_ptr()->exclusive_subjets_up_to(*this, nsub);
4024 }
4025 std::vector<PseudoJet> PseudoJet::exclusive_subjets (int nsub) const {
4026  vector<PseudoJet> subjets = exclusive_subjets_up_to(nsub);
4027  if (int(subjets.size()) < nsub) {
4028  ostringstream err;
4029  err << "Requested " << nsub << " exclusive subjets, but there were only "
4030  << subjets.size() << " particles in the jet";
4031  throw Error(err.str());
4032  }
4033  return subjets;
4034 }
4035 double PseudoJet::exclusive_subdmerge(int nsub) const {
4036  return validated_structure_ptr()->exclusive_subdmerge(*this, nsub);
4037 }
4038 double PseudoJet::exclusive_subdmerge_max(int nsub) const {
4039  return validated_structure_ptr()->exclusive_subdmerge_max(*this, nsub);
4040 }
4042  return ((_structure) && (_structure->has_pieces(*this)));
4043 }
4044 std::vector<PseudoJet> PseudoJet::pieces() const{
4045  return validated_structure_ptr()->pieces(*this);
4046 }
4047 PseudoJet::InexistentUserInfo::InexistentUserInfo() : Error("you attempted to perform a dynamic cast of a PseudoJet's extra info, but the extra info pointer was null")
4048 {}
4049 void sort_indices(vector<int> & indices,
4050  const vector<double> & values) {
4051  IndexedSortHelper index_sort_helper(&values);
4052  sort(indices.begin(), indices.end(), index_sort_helper);
4053 }
4054 vector<PseudoJet> sorted_by_pt(const vector<PseudoJet> & jets) {
4055  vector<double> minus_kt2(jets.size());
4056  for (size_t i = 0; i < jets.size(); i++) {minus_kt2[i] = -jets[i].kt2();}
4057  return objects_sorted_by_values(jets, minus_kt2);
4058 }
4059 vector<PseudoJet> sorted_by_rapidity(const vector<PseudoJet> & jets) {
4060  vector<double> rapidities(jets.size());
4061  for (size_t i = 0; i < jets.size(); i++) {rapidities[i] = jets[i].rap();}
4062  return objects_sorted_by_values(jets, rapidities);
4063 }
4064 vector<PseudoJet> sorted_by_E(const vector<PseudoJet> & jets) {
4065  vector<double> energies(jets.size());
4066  for (size_t i = 0; i < jets.size(); i++) {energies[i] = -jets[i].E();}
4067  return objects_sorted_by_values(jets, energies);
4068 }
4069 vector<PseudoJet> sorted_by_pz(const vector<PseudoJet> & jets) {
4070  vector<double> pz(jets.size());
4071  for (size_t i = 0; i < jets.size(); i++) {pz[i] = jets[i].pz();}
4072  return objects_sorted_by_values(jets, pz);
4073 }
4074 PseudoJet join(const vector<PseudoJet> & pieces){
4075  PseudoJet result; // automatically initialised to 0
4076  for (unsigned int i=0; i<pieces.size(); i++)
4077  result += pieces[i];
4078  CompositeJetStructure *cj_struct = new CompositeJetStructure(pieces);
4079  result.set_structure_shared_ptr(SharedPtr<PseudoJetStructureBase>(cj_struct));
4080  return result;
4081 }
4083  return join(vector<PseudoJet>(1,j1));
4084 }
4085 PseudoJet join(const PseudoJet & j1, const PseudoJet & j2){
4086  vector<PseudoJet> pieces;
4087  pieces.reserve(2);
4088  pieces.push_back(j1);
4089  pieces.push_back(j2);
4090  return join(pieces);
4091 }
4092 PseudoJet join(const PseudoJet & j1, const PseudoJet & j2, const PseudoJet & j3){
4093  vector<PseudoJet> pieces;
4094  pieces.reserve(3);
4095  pieces.push_back(j1);
4096  pieces.push_back(j2);
4097  pieces.push_back(j3);
4098  return join(pieces);
4099 }
4100 PseudoJet join(const PseudoJet & j1, const PseudoJet & j2, const PseudoJet & j3, const PseudoJet & j4){
4101  vector<PseudoJet> pieces;
4102  pieces.reserve(4);
4103  pieces.push_back(j1);
4104  pieces.push_back(j2);
4105  pieces.push_back(j3);
4106  pieces.push_back(j4);
4107  return join(pieces);
4108 }
4110 using namespace std;
4111 FJCORE_BEGIN_NAMESPACE // defined in fastjet/internal/base.hh
4113  return NULL;
4114 }
4115 const ClusterSequence * PseudoJetStructureBase::validated_cs() const{
4116  throw Error("This PseudoJet structure is not associated with a valid ClusterSequence");
4117 }
4118 bool PseudoJetStructureBase::has_partner(const PseudoJet & /*reference */, PseudoJet & /*partner*/) const{
4119  throw Error("This PseudoJet structure has no implementation for has_partner");
4120 }
4121 bool PseudoJetStructureBase::has_child(const PseudoJet & /*reference*/, PseudoJet & /*child*/) const{
4122  throw Error("This PseudoJet structure has no implementation for has_child");
4123 }
4124 bool PseudoJetStructureBase::has_parents(const PseudoJet & /*reference*/, PseudoJet &/*parent1*/, PseudoJet &/*parent2*/) const{
4125  throw Error("This PseudoJet structure has no implementation for has_parents");
4126 }
4127 bool PseudoJetStructureBase::object_in_jet(const PseudoJet & /*reference*/, const PseudoJet & /*jet*/) const{
4128  throw Error("This PseudoJet structure has no implementation for is_inside");
4129 }
4130 vector<PseudoJet> PseudoJetStructureBase::constituents(const PseudoJet &/*reference*/) const{
4131  throw Error("This PseudoJet structure has no implementation for constituents");
4132 }
4133 vector<PseudoJet> PseudoJetStructureBase::exclusive_subjets (const PseudoJet & /*reference*/, const double & /*dcut*/) const{
4134  throw Error("This PseudoJet structure has no implementation for exclusive_subjets");
4135 }
4136 int PseudoJetStructureBase::n_exclusive_subjets(const PseudoJet & /*reference*/, const double & /*dcut*/) const{
4137  throw Error("This PseudoJet structure has no implementation for n_exclusive_subjets");
4138 }
4139 vector<PseudoJet> PseudoJetStructureBase::exclusive_subjets_up_to (const PseudoJet & /*reference*/, int /*nsub*/) const{
4140  throw Error("This PseudoJet structure has no implementation for exclusive_subjets");
4141 }
4142 double PseudoJetStructureBase::exclusive_subdmerge(const PseudoJet & /*reference*/, int /*nsub*/) const{
4143  throw Error("This PseudoJet structure has no implementation for exclusive_submerge");
4144 }
4145 double PseudoJetStructureBase::exclusive_subdmerge_max(const PseudoJet & /*reference*/, int /*nsub*/) const{
4146  throw Error("This PseudoJet structure has no implementation for exclusive_submerge_max");
4147 }
4148 std::vector<PseudoJet> PseudoJetStructureBase::pieces(const PseudoJet & /*reference*/) const{
4149  throw Error("This PseudoJet structure has no implementation for pieces");
4150 }
4152 #include <sstream>
4153 #include <algorithm>
4154 using namespace std;
4155 FJCORE_BEGIN_NAMESPACE // defined in fastjet/internal/base.hh
4156 std::vector<PseudoJet> Selector::operator()(const std::vector<PseudoJet> & jets) const {
4157  std::vector<PseudoJet> result;
4158  const SelectorWorker * worker_local = validated_worker();
4159  if (worker_local->applies_jet_by_jet()) {
4160  for (std::vector<PseudoJet>::const_iterator jet = jets.begin();
4161  jet != jets.end(); jet++) {
4162  if (worker_local->pass(*jet)) result.push_back(*jet);
4163  }
4164  } else {
4165  std::vector<const PseudoJet *> jetptrs(jets.size());
4166  for (unsigned i = 0; i < jets.size(); i++) {
4167  jetptrs[i] = & jets[i];
4168  }
4169  worker_local->terminator(jetptrs);
4170  for (unsigned i = 0; i < jetptrs.size(); i++) {
4171  if (jetptrs[i]) result.push_back(jets[i]);
4172  }
4173  }
4174  return result;
4175 }
4176 unsigned int Selector::count(const std::vector<PseudoJet> & jets) const {
4177  unsigned n = 0;
4178  const SelectorWorker * worker_local = validated_worker();
4179  if (worker_local->applies_jet_by_jet()) {
4180  for (unsigned i = 0; i < jets.size(); i++) {
4181  if (worker_local->pass(jets[i])) n++;
4182  }
4183  } else {
4184  std::vector<const PseudoJet *> jetptrs(jets.size());
4185  for (unsigned i = 0; i < jets.size(); i++) {
4186  jetptrs[i] = & jets[i];
4187  }
4188  worker_local->terminator(jetptrs);
4189  for (unsigned i = 0; i < jetptrs.size(); i++) {
4190  if (jetptrs[i]) n++;
4191  }
4192  }
4193  return n;
4194 }
4195 PseudoJet Selector::sum(const std::vector<PseudoJet> & jets) const {
4196  PseudoJet this_sum(0,0,0,0);
4197  const SelectorWorker * worker_local = validated_worker();
4198  if (worker_local->applies_jet_by_jet()) {
4199  for (unsigned i = 0; i < jets.size(); i++) {
4200  if (worker_local->pass(jets[i])) this_sum += jets[i];
4201  }
4202  } else {
4203  std::vector<const PseudoJet *> jetptrs(jets.size());
4204  for (unsigned i = 0; i < jets.size(); i++) {
4205  jetptrs[i] = & jets[i];
4206  }
4207  worker_local->terminator(jetptrs);
4208  for (unsigned i = 0; i < jetptrs.size(); i++) {
4209  if (jetptrs[i]) this_sum += jets[i];
4210  }
4211  }
4212  return this_sum;
4213 }
4214 double Selector::scalar_pt_sum(const std::vector<PseudoJet> & jets) const {
4215  double this_sum = 0.0;
4216  const SelectorWorker * worker_local = validated_worker();
4217  if (worker_local->applies_jet_by_jet()) {
4218  for (unsigned i = 0; i < jets.size(); i++) {
4219  if (worker_local->pass(jets[i])) this_sum += jets[i].pt();
4220  }
4221  } else {
4222  std::vector<const PseudoJet *> jetptrs(jets.size());
4223  for (unsigned i = 0; i < jets.size(); i++) {
4224  jetptrs[i] = & jets[i];
4225  }
4226  worker_local->terminator(jetptrs);
4227  for (unsigned i = 0; i < jetptrs.size(); i++) {
4228  if (jetptrs[i]) this_sum += jets[i].pt();
4229  }
4230  }
4231  return this_sum;
4232 }
4233 void Selector::sift(const std::vector<PseudoJet> & jets,
4234  std::vector<PseudoJet> & jets_that_pass,
4235  std::vector<PseudoJet> & jets_that_fail
4236  ) const {
4237  const SelectorWorker * worker_local = validated_worker();
4238  jets_that_pass.clear();
4239  jets_that_fail.clear();
4240  if (worker_local->applies_jet_by_jet()) {
4241  for (unsigned i = 0; i < jets.size(); i++) {
4242  if (worker_local->pass(jets[i])) {
4243  jets_that_pass.push_back(jets[i]);
4244  } else {
4245  jets_that_fail.push_back(jets[i]);
4246  }
4247  }
4248  } else {
4249  std::vector<const PseudoJet *> jetptrs(jets.size());
4250  for (unsigned i = 0; i < jets.size(); i++) {
4251  jetptrs[i] = & jets[i];
4252  }
4253  worker_local->terminator(jetptrs);
4254  for (unsigned i = 0; i < jetptrs.size(); i++) {
4255  if (jetptrs[i]) {
4256  jets_that_pass.push_back(jets[i]);
4257  } else {
4258  jets_that_fail.push_back(jets[i]);
4259  }
4260  }
4261  }
4262 }
4264  if (! is_geometric()) return false;
4265  double rapmin, rapmax;
4266  get_rapidity_extent(rapmin, rapmax);
4267  return (rapmax != std::numeric_limits<double>::infinity())
4268  && (-rapmin != std::numeric_limits<double>::infinity());
4269 }
4270 class SW_Identity : public SelectorWorker {
4271 public:
4273  virtual bool pass(const PseudoJet &) const {
4274  return true;
4275  }
4276  virtual void terminator(vector<const PseudoJet *> &) const {
4277  return;
4278  }
4279  virtual string description() const { return "Identity";}
4280  virtual bool is_geometric() const { return true;}
4281 };
4283  return Selector(new SW_Identity);
4284 }
4285 class SW_Not : public SelectorWorker {
4286 public:
4287  SW_Not(const Selector & s) : _s(s) {}
4288  virtual SelectorWorker* copy(){ return new SW_Not(*this);}
4289  virtual bool pass(const PseudoJet & jet) const {
4290  if (!applies_jet_by_jet())
4291  throw Error("Cannot apply this selector worker to an individual jet");
4292  return ! _s.pass(jet);
4293  }
4294  virtual bool applies_jet_by_jet() const {return _s.applies_jet_by_jet();}
4295  virtual void terminator(vector<const PseudoJet *> & jets) const {
4296  if (applies_jet_by_jet()){
4298  return;
4299  }
4300  vector<const PseudoJet *> s_jets = jets;
4301  _s.worker()->terminator(s_jets);
4302  for (unsigned int i=0; i<s_jets.size(); i++){
4303  if (s_jets[i]) jets[i] = NULL;
4304  }
4305  }
4306  virtual string description() const {
4307  ostringstream ostr;
4308  ostr << "!(" << _s.description() << ")";
4309  return ostr.str();
4310  }
4311  virtual bool is_geometric() const { return _s.is_geometric();}
4312  virtual bool takes_reference() const { return _s.takes_reference();}
4313  virtual void set_reference(const PseudoJet &ref) { _s.set_reference(ref);}
4314 protected:
4316 };
4318  return Selector(new SW_Not(s));
4319 }
4321 public:
4322  SW_BinaryOperator(const Selector & s1, const Selector & s2) : _s1(s1), _s2(s2) {
4323  _applies_jet_by_jet = _s1.applies_jet_by_jet() && _s2.applies_jet_by_jet();
4324  _takes_reference = _s1.takes_reference() || _s2.takes_reference();
4325  _is_geometric = _s1.is_geometric() && _s2.is_geometric();
4326  }
4327  virtual bool applies_jet_by_jet() const {return _applies_jet_by_jet;}
4328  virtual bool takes_reference() const{
4329  return _takes_reference;
4330  }
4331  virtual void set_reference(const PseudoJet &centre){
4332  _s1.set_reference(centre);
4333  _s2.set_reference(centre);
4334  }
4335  virtual bool is_geometric() const { return _is_geometric;}
4336 protected:
4341 };
4342 class SW_And: public SW_BinaryOperator {
4343 public:
4344  SW_And(const Selector & s1, const Selector & s2) : SW_BinaryOperator(s1,s2){}
4345  virtual SelectorWorker* copy(){ return new SW_And(*this);}
4346  virtual bool pass(const PseudoJet & jet) const {
4347  if (!applies_jet_by_jet())
4348  throw Error("Cannot apply this selector worker to an individual jet");
4349  return _s1.pass(jet) && _s2.pass(jet);
4350  }
4351  virtual void terminator(vector<const PseudoJet *> & jets) const {
4352  if (applies_jet_by_jet()){
4354  return;
4355  }
4356  vector<const PseudoJet *> s1_jets = jets;
4357  _s1.worker()->terminator(s1_jets);
4358  _s2.worker()->terminator(jets);
4359  for (unsigned int i=0; i<jets.size(); i++){
4360  if (! s1_jets[i]) jets[i] = NULL;
4361  }
4362  }
4363  virtual void get_rapidity_extent(double & rapmin, double & rapmax) const {
4364  double s1min, s1max, s2min, s2max;
4365  _s1.get_rapidity_extent(s1min, s1max);
4366  _s2.get_rapidity_extent(s2min, s2max);
4367  rapmax = min(s1max, s2max);
4368  rapmin = max(s1min, s2min);
4369  }
4370  virtual string description() const {
4371  ostringstream ostr;
4372  ostr << "(" << _s1.description() << " && " << _s2.description() << ")";
4373  return ostr.str();
4374  }
4375 };
4376 Selector operator&&(const Selector & s1, const Selector & s2) {
4377  return Selector(new SW_And(s1,s2));
4378 }
4379 class SW_Or: public SW_BinaryOperator {
4380 public:
4381  SW_Or(const Selector & s1, const Selector & s2) : SW_BinaryOperator(s1,s2) {}
4382  virtual SelectorWorker* copy(){ return new SW_Or(*this);}
4383  virtual bool pass(const PseudoJet & jet) const {
4384  if (!applies_jet_by_jet())
4385  throw Error("Cannot apply this selector worker to an individual jet");
4386  return _s1.pass(jet) || _s2.pass(jet);
4387  }
4388  virtual bool applies_jet_by_jet() const {
4389  return _s1.applies_jet_by_jet() && _s2.applies_jet_by_jet();
4390  }
4391  virtual void terminator(vector<const PseudoJet *> & jets) const {
4392  if (applies_jet_by_jet()){
4394  return;
4395  }
4396  vector<const PseudoJet *> s1_jets = jets;
4397  _s1.worker()->terminator(s1_jets);
4398  _s2.worker()->terminator(jets);
4399  for (unsigned int i=0; i<jets.size(); i++){
4400  if (s1_jets[i]) jets[i] = s1_jets[i];
4401  }
4402  }
4403  virtual string description() const {
4404  ostringstream ostr;
4405  ostr << "(" << _s1.description() << " || " << _s2.description() << ")";
4406  return ostr.str();
4407  }
4408  virtual void get_rapidity_extent(double & rapmin, double & rapmax) const {
4409  double s1min, s1max, s2min, s2max;
4410  _s1.get_rapidity_extent(s1min, s1max);
4411  _s2.get_rapidity_extent(s2min, s2max);
4412  rapmax = max(s1max, s2max);
4413  rapmin = min(s1min, s2min);
4414  }
4415 };
4416 Selector operator ||(const Selector & s1, const Selector & s2) {
4417  return Selector(new SW_Or(s1,s2));
4418 }
4419 class SW_Mult: public SW_And {
4420 public:
4421  SW_Mult(const Selector & s1, const Selector & s2) : SW_And(s1,s2) {}
4422  virtual SelectorWorker* copy(){ return new SW_Mult(*this);}
4423  virtual void terminator(vector<const PseudoJet *> & jets) const {
4424  if (applies_jet_by_jet()){
4426  return;
4427  }
4428  _s2.worker()->terminator(jets);
4429  _s1.worker()->terminator(jets);
4430  }
4431  virtual string description() const {
4432  ostringstream ostr;
4433  ostr << "(" << _s1.description() << " * " << _s2.description() << ")";
4434  return ostr.str();
4435  }
4436 };
4437 Selector operator*(const Selector & s1, const Selector & s2) {
4438  return Selector(new SW_Mult(s1,s2));
4439 }
4441 public:
4442  QuantityBase(double q) : _q(q){}
4443  virtual ~QuantityBase(){}
4444  virtual double operator()(const PseudoJet & jet ) const =0;
4445  virtual string description() const =0;
4446  virtual bool is_geometric() const { return false;}
4447  virtual double comparison_value() const {return _q;}
4448  virtual double description_value() const {return comparison_value();}
4449 protected:
4450  double _q;
4451 };
4453 public:
4454  QuantitySquareBase(double sqrtq) : QuantityBase(sqrtq*sqrtq), _sqrtq(sqrtq){}
4455  virtual double description_value() const {return _sqrtq;}
4456 protected:
4457  double _sqrtq;
4458 };
4459 template<typename QuantityType>
4460 class SW_QuantityMin : public SelectorWorker{
4461 public:
4462  SW_QuantityMin(double qmin) : _qmin(qmin) {}
4463  virtual bool pass(const PseudoJet & jet) const {return _qmin(jet) >= _qmin.comparison_value();}
4464  virtual string description() const {
4465  ostringstream ostr;
4466  ostr << _qmin.description() << " >= " << _qmin.description_value();
4467  return ostr.str();
4468  }
4469  virtual bool is_geometric() const { return _qmin.is_geometric();}
4470 protected:
4471  QuantityType _qmin;
4472 };
4473 template<typename QuantityType>
4474 class SW_QuantityMax : public SelectorWorker {
4475 public:
4476  SW_QuantityMax(double qmax) : _qmax(qmax) {}
4477  virtual bool pass(const PseudoJet & jet) const {return _qmax(jet) <= _qmax.comparison_value();}
4478  virtual string description() const {
4479  ostringstream ostr;
4480  ostr << _qmax.description() << " <= " << _qmax.description_value();
4481  return ostr.str();
4482  }
4483  virtual bool is_geometric() const { return _qmax.is_geometric();}
4484 protected:
4485  QuantityType _qmax;
4486 };
4487 template<typename QuantityType>
4489 public:
4490  SW_QuantityRange(double qmin, double qmax) : _qmin(qmin), _qmax(qmax) {}
4491  virtual bool pass(const PseudoJet & jet) const {
4492  double q = _qmin(jet); // we could identically use _qmax
4493  return (q >= _qmin.comparison_value()) && (q <= _qmax.comparison_value());
4494  }
4495  virtual string description() const {
4496  ostringstream ostr;
4497  ostr << _qmin.description_value() << " <= " << _qmin.description() << " <= " << _qmax.description_value();
4498  return ostr.str();
4499  }
4500  virtual bool is_geometric() const { return _qmin.is_geometric();}
4501 protected:
4502  QuantityType _qmin; // the lower cut
4503  QuantityType _qmax; // the upper cut
4504 };
4506 public:
4508  virtual double operator()(const PseudoJet & jet ) const { return jet.perp2();}
4509  virtual string description() const {return "pt";}
4510 };
4511 Selector SelectorPtMin(double ptmin) {
4512  return Selector(new SW_QuantityMin<QuantityPt2>(ptmin));
4513 }
4514 Selector SelectorPtMax(double ptmax) {
4515  return Selector(new SW_QuantityMax<QuantityPt2>(ptmax));
4516 }
4517 Selector SelectorPtRange(double ptmin, double ptmax) {
4518  return Selector(new SW_QuantityRange<QuantityPt2>(ptmin, ptmax));
4519 }
4521 public:
4523  virtual double operator()(const PseudoJet & jet ) const { return jet.Et2();}
4524  virtual string description() const {return "Et";}
4525 };
4526 Selector SelectorEtMin(double Etmin) {
4527  return Selector(new SW_QuantityMin<QuantityEt2>(Etmin));
4528 }
4529 Selector SelectorEtMax(double Etmax) {
4530  return Selector(new SW_QuantityMax<QuantityEt2>(Etmax));
4531 }
4532 Selector SelectorEtRange(double Etmin, double Etmax) {
4533  return Selector(new SW_QuantityRange<QuantityEt2>(Etmin, Etmax));
4534 }
4535 class QuantityE : public QuantityBase{
4536 public:
4537  QuantityE(double E) : QuantityBase(E){}
4538  virtual double operator()(const PseudoJet & jet ) const { return jet.E();}
4539  virtual string description() const {return "E";}
4540 };
4541 Selector SelectorEMin(double Emin) {
4542  return Selector(new SW_QuantityMin<QuantityE>(Emin));
4543 }
4544 Selector SelectorEMax(double Emax) {
4545  return Selector(new SW_QuantityMax<QuantityE>(Emax));
4546 }
4547 Selector SelectorERange(double Emin, double Emax) {
4548  return Selector(new SW_QuantityRange<QuantityE>(Emin, Emax));
4549 }
4551 public:
4553  virtual double operator()(const PseudoJet & jet ) const { return jet.m2();}
4554  virtual string description() const {return "mass";}
4555 };
4557  return Selector(new SW_QuantityMin<QuantityM2>(mmin));
4558 }
4560  return Selector(new SW_QuantityMax<QuantityM2>(mmax));
4561 }
4562 Selector SelectorMassRange(double mmin, double mmax) {
4563  return Selector(new SW_QuantityRange<QuantityM2>(mmin, mmax));
4564 }
4565 class QuantityRap : public QuantityBase{
4566 public:
4567  QuantityRap(double rap) : QuantityBase(rap){}
4568  virtual double operator()(const PseudoJet & jet ) const { return jet.rap();}
4569  virtual string description() const {return "rap";}
4570  virtual bool is_geometric() const { return true;}
4571 };
4572 class SW_RapMin : public SW_QuantityMin<QuantityRap>{
4573 public:
4574  SW_RapMin(double rapmin) : SW_QuantityMin<QuantityRap>(rapmin){}
4575  virtual void get_rapidity_extent(double &rapmin, double & rapmax) const{
4576  rapmax = std::numeric_limits<double>::max();
4577  rapmin = _qmin.comparison_value();
4578  }
4579 };
4580 class SW_RapMax : public SW_QuantityMax<QuantityRap>{
4581 public:
4582  SW_RapMax(double rapmax) : SW_QuantityMax<QuantityRap>(rapmax){}
4583  virtual void get_rapidity_extent(double &rapmin, double & rapmax) const{
4584  rapmax = _qmax.comparison_value();
4585  rapmin = -std::numeric_limits<double>::max();
4586  }
4587 };
4588 class SW_RapRange : public SW_QuantityRange<QuantityRap>{
4589 public:
4590  SW_RapRange(double rapmin, double rapmax) : SW_QuantityRange<QuantityRap>(rapmin, rapmax){
4591  assert(rapmin<=rapmax);
4592  }
4593  virtual void get_rapidity_extent(double &rapmin, double & rapmax) const{
4594  rapmax = _qmax.comparison_value();
4595  rapmin = _qmin.comparison_value();
4596  }
4597  virtual bool has_known_area() const { return true;}
4598  virtual double known_area() const {
4599  return twopi * (_qmax.comparison_value()-_qmin.comparison_value());
4600  }
4601 };
4602 Selector SelectorRapMin(double rapmin) {
4603  return Selector(new SW_RapMin(rapmin));
4604 }
4605 Selector SelectorRapMax(double rapmax) {
4606  return Selector(new SW_RapMax(rapmax));
4607 }
4608 Selector SelectorRapRange(double rapmin, double rapmax) {
4609  return Selector(new SW_RapRange(rapmin, rapmax));
4610 }
4612 public:
4613  QuantityAbsRap(double absrap) : QuantityBase(absrap){}
4614  virtual double operator()(const PseudoJet & jet ) const { return abs(jet.rap());}
4615  virtual string description() const {return "|rap|";}
4616  virtual bool is_geometric() const { return true;}
4617 };
4618 class SW_AbsRapMax : public SW_QuantityMax<QuantityAbsRap>{
4619 public:
4620  SW_AbsRapMax(double absrapmax) : SW_QuantityMax<QuantityAbsRap>(absrapmax){}
4621  virtual void get_rapidity_extent(double &rapmin, double & rapmax) const{
4622  rapmax = _qmax.comparison_value();
4623  rapmin = -_qmax.comparison_value();
4624  }
4625  virtual bool has_known_area() const { return true;}
4626  virtual double known_area() const {
4627  return twopi * 2 * _qmax.comparison_value();
4628  }
4629 };
4630 class SW_AbsRapRange : public SW_QuantityRange<QuantityAbsRap>{
4631 public:
4632  SW_AbsRapRange(double absrapmin, double absrapmax) : SW_QuantityRange<QuantityAbsRap>(absrapmin, absrapmax){}
4633  virtual void get_rapidity_extent(double &rapmin, double & rapmax) const{
4634  rapmax = _qmax.comparison_value();
4635  rapmin = -_qmax.comparison_value();
4636  }
4637  virtual bool has_known_area() const { return true;}
4638  virtual double known_area() const {
4639  return twopi * 2 * (_qmax.comparison_value()-max(_qmin.comparison_value(),0.0)); // this should handle properly absrapmin<0
4640  }
4641 };
4642 Selector SelectorAbsRapMin(double absrapmin) {
4643  return Selector(new SW_QuantityMin<QuantityAbsRap>(absrapmin));
4644 }
4645 Selector SelectorAbsRapMax(double absrapmax) {
4646  return Selector(new SW_AbsRapMax(absrapmax));
4647 }
4648 Selector SelectorAbsRapRange(double rapmin, double rapmax) {
4649  return Selector(new SW_AbsRapRange(rapmin, rapmax));
4650 }
4651 class QuantityEta : public QuantityBase{
4652 public:
4653  QuantityEta(double eta) : QuantityBase(eta){}
4654  virtual double operator()(const PseudoJet & jet ) const { return jet.eta();}
4655  virtual string description() const {return "eta";}
4656 };
4658  return Selector(new SW_QuantityMin<QuantityEta>(etamin));
4659 }
4661  return Selector(new SW_QuantityMax<QuantityEta>(etamax));
4662 }
4664  return Selector(new SW_QuantityRange<QuantityEta>(etamin, etamax));
4665 }
4667 public:
4668  QuantityAbsEta(double abseta) : QuantityBase(abseta){}
4669  virtual double operator()(const PseudoJet & jet ) const { return abs(jet.eta());}
4670  virtual string description() const {return "|eta|";}
4671  virtual bool is_geometric() const { return true;}
4672 };
4673 Selector SelectorAbsEtaMin(double absetamin) {
4674  return Selector(new SW_QuantityMin<QuantityAbsEta>(absetamin));
4675 }
4676 Selector SelectorAbsEtaMax(double absetamax) {
4677  return Selector(new SW_QuantityMax<QuantityAbsEta>(absetamax));
4678 }
4679 Selector SelectorAbsEtaRange(double absetamin, double absetamax) {
4680  return Selector(new SW_QuantityRange<QuantityAbsEta>(absetamin, absetamax));
4681 }
4682 class SW_PhiRange : public SelectorWorker {
4683 public:
4684  SW_PhiRange(double phimin, double phimax) : _phimin(phimin), _phimax(phimax){
4685  assert(_phimin<_phimax);
4686  assert(_phimin>-twopi);
4687  assert(_phimax<2*twopi);
4688  _phispan = _phimax - _phimin;
4689  }
4690  virtual bool pass(const PseudoJet & jet) const {
4691  double dphi=jet.phi()-_phimin;
4692  if (dphi >= twopi) dphi -= twopi;
4693  if (dphi < 0) dphi += twopi;
4694  return (dphi <= _phispan);
4695  }
4696  virtual string description() const {
4697  ostringstream ostr;
4698  ostr << _phimin << " <= phi <= " << _phimax;
4699  return ostr.str();
4700  }
4701  virtual bool is_geometric() const { return true;}
4702 protected:
4703  double _phimin; // the lower cut
4704  double _phimax; // the upper cut
4705  double _phispan; // the span of the range
4706 };
4707 Selector SelectorPhiRange(double phimin, double phimax) {
4708  return Selector(new SW_PhiRange(phimin, phimax));
4709 }
4710 class SW_RapPhiRange : public SW_And{
4711 public:
4712  SW_RapPhiRange(double rapmin, double rapmax, double phimin, double phimax)
4713  : SW_And(SelectorRapRange(rapmin, rapmax), SelectorPhiRange(phimin, phimax)){
4714  _known_area = ((phimax-phimin > twopi) ? twopi : phimax-phimin) * (rapmax-rapmin);
4715  }
4716  virtual double known_area() const{
4717  return _known_area;
4718  }
4719 protected:
4720  double _known_area;
4721 };
4722 Selector SelectorRapPhiRange(double rapmin, double rapmax, double phimin, double phimax) {
4723  return Selector(new SW_RapPhiRange(rapmin, rapmax, phimin, phimax));
4724 }
4725 class SW_NHardest : public SelectorWorker {
4726 public:
4727  SW_NHardest(unsigned int n) : _n(n) {};
4728  virtual bool pass(const PseudoJet &) const {
4729  if (!applies_jet_by_jet())
4730  throw Error("Cannot apply this selector worker to an individual jet");
4731  return false;
4732  }
4733  virtual void terminator(vector<const PseudoJet *> & jets) const {
4734  if (jets.size() < _n) return;
4735  vector<double> minus_pt2(jets.size());
4736  vector<unsigned int> indices(jets.size());
4737  for (unsigned int i=0; i<jets.size(); i++){
4738  indices[i] = i;
4739  minus_pt2[i] = jets[i] ? -jets[i]->perp2() : 0.0;
4740  }
4741  IndexedSortHelper sort_helper(& minus_pt2);
4742  partial_sort(indices.begin(), indices.begin()+_n, indices.end(), sort_helper);
4743  for (unsigned int i=_n; i<jets.size(); i++)
4744  jets[indices[i]] = NULL;
4745  }
4746  virtual bool applies_jet_by_jet() const {return false;}
4747  virtual string description() const {
4748  ostringstream ostr;
4749  ostr << _n << " hardest";
4750  return ostr.str();
4751  }
4752 protected:
4753  unsigned int _n;
4754 };
4755 Selector SelectorNHardest(unsigned int n) {
4756  return Selector(new SW_NHardest(n));
4757 }
4758 class SW_WithReference : public SelectorWorker{
4759 public:
4760  SW_WithReference() : _is_initialised(false){};
4761  virtual bool takes_reference() const { return true;}
4762  virtual void set_reference(const PseudoJet &centre){
4763  _is_initialised = true;
4764  _reference = centre;
4765  }
4766 protected:
4769 };
4770 class SW_Circle : public SW_WithReference {
4771 public:
4772  SW_Circle(const double radius) : _radius2(radius*radius) {}
4773  virtual SelectorWorker* copy(){ return new SW_Circle(*this);}
4774  virtual bool pass(const PseudoJet & jet) const {
4775  if (! _is_initialised)
4776  throw Error("To use a SelectorCircle (or any selector that requires a reference), you first have to call set_reference(...)");
4777  return jet.squared_distance(_reference) <= _radius2;
4778  }
4779  virtual string description() const {
4780  ostringstream ostr;
4781  ostr << "distance from the centre <= " << sqrt(_radius2);
4782  return ostr.str();
4783  }
4784  virtual void get_rapidity_extent(double & rapmin, double & rapmax) const{
4785  if (! _is_initialised)
4786  throw Error("To use a SelectorCircle (or any selector that requires a reference), you first have to call set_reference(...)");
4787  rapmax = _reference.rap()+sqrt(_radius2);
4788  rapmin = _reference.rap()-sqrt(_radius2);
4789  }
4790  virtual bool is_geometric() const { return true;}
4791  virtual bool has_finite_area() const { return true;}
4792  virtual bool has_known_area() const { return true;}
4793  virtual double known_area() const {
4794  return pi * _radius2;
4795  }
4796 protected:
4797  double _radius2;
4798 };
4799 Selector SelectorCircle(const double radius) {
4800  return Selector(new SW_Circle(radius));
4801 }
4803 public:
4804  SW_Doughnut(const double radius_in, const double radius_out)
4805  : _radius_in2(radius_in*radius_in), _radius_out2(radius_out*radius_out) {}
4806  virtual SelectorWorker* copy(){ return new SW_Doughnut(*this);}
4807  virtual bool pass(const PseudoJet & jet) const {
4808  if (! _is_initialised)
4809  throw Error("To use a SelectorDoughnut (or any selector that requires a reference), you first have to call set_reference(...)");
4810  double distance2 = jet.squared_distance(_reference);
4811  return (distance2 <= _radius_out2) && (distance2 >= _radius_in2);
4812  }
4813  virtual string description() const {
4814  ostringstream ostr;
4815  ostr << sqrt(_radius_in2) << " <= distance from the centre <= " << sqrt(_radius_out2);
4816  return ostr.str();
4817  }
4818  virtual void get_rapidity_extent(double & rapmin, double & rapmax) const{
4819  if (! _is_initialised)
4820  throw Error("To use a SelectorDoughnut (or any selector that requires a reference), you first have to call set_reference(...)");
4821  rapmax = _reference.rap()+sqrt(_radius_out2);
4822  rapmin = _reference.rap()-sqrt(_radius_out2);
4823  }
4824  virtual bool is_geometric() const { return true;}
4825  virtual bool has_finite_area() const { return true;}
4826  virtual bool has_known_area() const { return true;}
4827  virtual double known_area() const {
4828  return pi * (_radius_out2-_radius_in2);
4829  }
4830 protected:
4831  double _radius_in2, _radius_out2;
4832 };
4833 Selector SelectorDoughnut(const double radius_in, const double radius_out) {
4834  return Selector(new SW_Doughnut(radius_in, radius_out));
4835 }
4836 class SW_Strip : public SW_WithReference {
4837 public:
4838  SW_Strip(const double delta) : _delta(delta) {}
4839  virtual SelectorWorker* copy(){ return new SW_Strip(*this);}
4840  virtual bool pass(const PseudoJet & jet) const {
4841  if (! _is_initialised)
4842  throw Error("To use a SelectorStrip (or any selector that requires a reference), you first have to call set_reference(...)");
4843  return abs(jet.rap()-_reference.rap()) <= _delta;
4844  }
4845  virtual string description() const {
4846  ostringstream ostr;
4847  ostr << "|rap - rap_reference| <= " << _delta;
4848  return ostr.str();
4849  }
4850  virtual void get_rapidity_extent(double & rapmin, double & rapmax) const{
4851  if (! _is_initialised)
4852  throw Error("To use a SelectorStrip (or any selector that requires a reference), you first have to call set_reference(...)");
4853  rapmax = _reference.rap()+_delta;
4854  rapmin = _reference.rap()-_delta;
4855  }
4856  virtual bool is_geometric() const { return true;}
4857  virtual bool has_finite_area() const { return true;}
4858  virtual bool has_known_area() const { return true;}
4859  virtual double known_area() const {
4860  return twopi * 2 * _delta;
4861  }
4862 protected:
4863  double _delta;
4864 };
4865 Selector SelectorStrip(const double half_width) {
4866  return Selector(new SW_Strip(half_width));
4867 }
4869 public:
4870  SW_Rectangle(const double delta_rap, const double delta_phi)
4871  : _delta_rap(delta_rap), _delta_phi(delta_phi) {}
4872  virtual SelectorWorker* copy(){ return new SW_Rectangle(*this);}
4873  virtual bool pass(const PseudoJet & jet) const {
4874  if (! _is_initialised)
4875  throw Error("To use a SelectorRectangle (or any selector that requires a reference), you first have to call set_reference(...)");
4876  return (abs(jet.rap()-_reference.rap()) <= _delta_rap) && (abs(jet.delta_phi_to(_reference)) <= _delta_phi);
4877  }
4878  virtual string description() const {
4879  ostringstream ostr;
4880  ostr << "|rap - rap_reference| <= " << _delta_rap << " && |phi - phi_reference| <= " << _delta_phi ;
4881  return ostr.str();
4882  }
4883  virtual void get_rapidity_extent(double & rapmin, double & rapmax) const{
4884  if (! _is_initialised)
4885  throw Error("To use a SelectorRectangle (or any selector that requires a reference), you first have to call set_reference(...)");
4886  rapmax = _reference.rap()+_delta_rap;
4887  rapmin = _reference.rap()-_delta_rap;
4888  }
4889  virtual bool is_geometric() const { return true;}
4890  virtual bool has_finite_area() const { return true;}
4891  virtual bool has_known_area() const { return true;}
4892  virtual double known_area() const {
4893  return 4 * _delta_rap * _delta_phi;
4894  }
4895 protected:
4896  double _delta_rap, _delta_phi;
4897 };
4898 Selector SelectorRectangle(const double half_rap_width, const double half_phi_width) {
4899  return Selector(new SW_Rectangle(half_rap_width, half_phi_width));
4900 }
4902 public:
4903  SW_PtFractionMin(double fraction) : _fraction2(fraction*fraction){}
4904  virtual SelectorWorker* copy(){ return new SW_PtFractionMin(*this);}
4905  virtual bool pass(const PseudoJet & jet) const {
4906  if (! _is_initialised)
4907  throw Error("To use a SelectorPtFractionMin (or any selector that requires a reference), you first have to call set_reference(...)");
4908  return (jet.perp2() >= _fraction2*_reference.perp2());
4909  }
4910  virtual string description() const {
4911  ostringstream ostr;
4912  ostr << "pt >= " << sqrt(_fraction2) << "* pt_ref";
4913  return ostr.str();
4914  }
4915 protected:
4916  double _fraction2;
4917 };
4919  return Selector(new SW_PtFractionMin(fraction));
4920 }
4921 class SW_IsZero : public SelectorWorker {
4922 public:
4924  virtual bool pass(const PseudoJet & jet) const {
4925  return jet==0;
4926  }
4927  virtual string description() const { return "zero";}
4928 };
4930  return Selector(new SW_IsZero());
4931 }
4933  _worker.reset(new SW_And(*this, b));
4934  return *this;
4935 }
4937  _worker.reset(new SW_Or(*this, b));
4938  return *this;
4939 }
4940 FJCORE_END_NAMESPACE // defined in fastjet/internal/base.hh
4941 #include <iomanip>
4942 using namespace std;
4943 FJCORE_BEGIN_NAMESPACE // defined in fastjet/internal/base.hh
4944 LazyTiling25::LazyTiling25(ClusterSequence & cs) :
4945  _cs(cs), _jets(cs.jets())
4946 {
4947 #ifdef INSTRUMENT2
4948  _ncall = 0; // gps tmp
4949  _ncall_dtt = 0; // gps tmp
4950 #endif // INSTRUMENT2
4951  _Rparam = cs.jet_def().R();
4952  _R2 = _Rparam * _Rparam;
4953  _invR2 = 1.0 / _R2;
4955 }
4957  double default_size = max(0.1,_Rparam)/2;
4958  _tile_size_eta = default_size;
4959  _n_tiles_phi = max(5,int(floor(twopi/default_size)));
4960  _tile_size_phi = twopi / _n_tiles_phi; // >= _Rparam and fits in 2pi
4961 #define _FJCORE_TILING25_USE_TILING_ANALYSIS_
4962 #ifdef _FASTJET_TILING25_USE_TILING_ANALYSIS_
4963  TilingExtent tiling_analysis(_cs);
4964  _tiles_eta_min = tiling_analysis.minrap();
4965  _tiles_eta_max = tiling_analysis.maxrap();
4966 #else // not _FASTJET_TILING25_USE_TILING_ANALYSIS_
4967  _tiles_eta_min = 0.0;
4968  _tiles_eta_max = 0.0;
4969  const double maxrap = 7.0;
4970  for(unsigned int i = 0; i < _jets.size(); i++) {
4971  double eta = _jets[i].rap();
4972  if (abs(eta) < maxrap) {
4973  if (eta < _tiles_eta_min) {_tiles_eta_min = eta;}
4974  if (eta > _tiles_eta_max) {_tiles_eta_max = eta;}
4975  }
4976  }
4977 #endif // _FASTJET_TILING25_USE_TILING_ANALYSIS_
4978 # define FJCORE_LAZY25_MIN3TILESY
4979 #ifdef FJCORE_LAZY25_MIN3TILESY
4982  _tiles_ieta_min = 0;
4983  _tiles_ieta_max = 2;
4985  } else {
4986 #endif //FASTJET_LAZY25_MIN3TILESY
4991 #ifdef FJCORE_LAZY25_MIN3TILESY
4992  }
4993 #endif
4996  vector<bool> use_periodic_delta_phi(_n_tiles_phi, false);
4997  if (_n_tiles_phi <= 5) {
4998  fill(use_periodic_delta_phi.begin(), use_periodic_delta_phi.end(), true);
4999  } else {
5000  use_periodic_delta_phi[0] = true;
5001  use_periodic_delta_phi[1] = true;
5002  use_periodic_delta_phi[_n_tiles_phi-2] = true;
5003  use_periodic_delta_phi[_n_tiles_phi-1] = true;
5004  }
5006  for (int ieta = _tiles_ieta_min; ieta <= _tiles_ieta_max; ieta++) {
5007  for (int iphi = 0; iphi < _n_tiles_phi; iphi++) {
5008  Tile25 * tile = & _tiles[_tile_index(ieta,iphi)];
5009  tile->head = NULL; // first element of tiles points to itself
5010  tile->begin_tiles[0] = tile;
5011  Tile25 ** pptile = & (tile->begin_tiles[0]);
5012  pptile++;
5013  tile->surrounding_tiles = pptile;
5014  if (ieta > _tiles_ieta_min) {
5015  // with the itile subroutine, we can safely run tiles from
5016  // idphi=-1 to idphi=+1, because it takes care of
5017  // negative and positive boundaries
5018  for (int idphi = -2; idphi <=+2; idphi++) {
5019  *pptile = & _tiles[_tile_index(ieta-1,iphi+idphi)];
5020  pptile++;
5021  }
5022  }
5023  if (ieta > _tiles_ieta_min + 1) {
5024  // with the itile subroutine, we can safely run tiles from
5025  // idphi=-1 to idphi=+1, because it takes care of
5026  // negative and positive boundaries
5027  for (int idphi = -2; idphi <= +2; idphi++) {
5028  *pptile = & _tiles[_tile_index(ieta-2,iphi+idphi)];
5029  pptile++;
5030  }
5031  }
5032  *pptile = & _tiles[_tile_index(ieta,iphi-1)];
5033  pptile++;
5034  *pptile = & _tiles[_tile_index(ieta,iphi-2)];
5035  pptile++;
5036  tile->RH_tiles = pptile;
5037  *pptile = & _tiles[_tile_index(ieta,iphi+1)];
5038  pptile++;
5039  *pptile = & _tiles[_tile_index(ieta,iphi+2)];
5040  pptile++;
5041  if (ieta < _tiles_ieta_max) {
5042  for (int idphi = -2; idphi <= +2; idphi++) {
5043  *pptile = & _tiles[_tile_index(ieta+1,iphi+idphi)];
5044  pptile++;
5045  }
5046  }
5047  if (ieta < _tiles_ieta_max - 1) {
5048  for (int idphi = -2; idphi <= +2; idphi++) {
5049  *pptile = & _tiles[_tile_index(ieta+2,iphi+idphi)];
5050  pptile++;
5051  }
5052  }
5053  tile->end_tiles = pptile;
5054  tile->tagged = false;
5055  tile->use_periodic_delta_phi = use_periodic_delta_phi[iphi];
5056  tile->max_NN_dist = 0;
5057  tile->eta_centre = (ieta-_tiles_ieta_min+0.5)*_tile_size_eta + _tiles_eta_min;
5058  tile->phi_centre = (iphi+0.5)*_tile_size_phi;
5059  }
5060  }
5061 }
5062 int LazyTiling25::_tile_index(const double eta, const double phi) const {
5063  int ieta, iphi;
5064  if (eta <= _tiles_eta_min) {ieta = 0;}
5065  else if (eta >= _tiles_eta_max) {ieta = _tiles_ieta_max-_tiles_ieta_min;}
5066  else {
5067  ieta = int(((eta - _tiles_eta_min) / _tile_size_eta));
5068  if (ieta > _tiles_ieta_max-_tiles_ieta_min) {
5070  }
5071  iphi = int((phi+twopi)/_tile_size_phi) % _n_tiles_phi;
5072  return (iphi + ieta * _n_tiles_phi);
5073 }
5074 inline void LazyTiling25::_tj_set_jetinfo( TiledJet * const jet,
5075  const int _jets_index) {
5076  _bj_set_jetinfo<>(jet, _jets_index);
5077  jet->tile_index = _tile_index(jet->eta, jet->phi);
5078  Tile25 * tile = &_tiles[jet->tile_index];
5079  jet->previous = NULL;
5080  jet->next = tile->head;
5081  if (jet->next != NULL) {jet->next->previous = jet;}
5082  tile->head = jet;
5083 }
5085  Tile25 * tile = & _tiles[jet->tile_index];
5086  if (jet->previous == NULL) {
5087  tile->head = jet->next;
5088  } else {
5089  jet->previous->next = jet->next;
5090  }
5091  if (jet->next != NULL) {
5092  jet->next->previous = jet->previous;
5093  }
5094 }
5095 void LazyTiling25::_print_tiles(TiledJet * briefjets ) const {
5096  for (vector<Tile25>::const_iterator tile = _tiles.begin();
5097  tile < _tiles.end(); tile++) {
5098  cout << "Tile " << tile - _tiles.begin()
5099  << " at " << setw(10) << tile->eta_centre << "," << setw(10) << tile->phi_centre
5100  << " = ";
5101  vector<int> list;
5102  for (TiledJet * jetI = tile->head; jetI != NULL; jetI = jetI->next) {
5103  list.push_back(jetI-briefjets);
5104  }
5105  sort(list.begin(),list.end());
5106  for (unsigned int i = 0; i < list.size(); i++) {cout <<" "<<list[i];}
5107  cout <<"\n";
5108  }
5109 }
5111  vector<int> & tile_union, int & n_near_tiles) const {
5112  for (Tile25 * const * near_tile = _tiles[tile_index].begin_tiles;
5113  near_tile != _tiles[tile_index].end_tiles; near_tile++){
5114  tile_union[n_near_tiles] = *near_tile - & _tiles[0];
5115  n_near_tiles++;
5116  }
5117 }
5119  const int tile_index,
5120  vector<int> & tile_union, int & n_near_tiles) {
5121  for (Tile25 ** near_tile = _tiles[tile_index].begin_tiles;
5122  near_tile != _tiles[tile_index].end_tiles; near_tile++){
5123  if (! (*near_tile)->tagged) {
5124  (*near_tile)->tagged = true;
5125  tile_union[n_near_tiles] = *near_tile - & _tiles[0];
5126  n_near_tiles++;
5127  }
5128  }
5129 }
5131  const TiledJet * jet,
5132  vector<int> & tile_union, int & n_near_tiles) {
5133  Tile25 & tile = _tiles[jet->tile_index];
5134  for (Tile25 ** near_tile = tile.begin_tiles; near_tile != tile.end_tiles; near_tile++){
5135  if ((*near_tile)->tagged) continue;
5136  double dist = _distance_to_tile(jet, *near_tile) - tile_edge_security_margin;
5137  if (dist > (*near_tile)->max_NN_dist) continue;
5138  (*near_tile)->tagged = true;
5139  tile_union[n_near_tiles] = *near_tile - & _tiles[0];
5140  n_near_tiles++;
5141  }
5142 }
5143 inline double LazyTiling25::_distance_to_tile(const TiledJet * bj, const Tile25 * tile)
5144 #ifdef INSTRUMENT2
5145  {
5146  _ncall_dtt++; // GPS tmp
5147 #else
5148  const {
5149 #endif // INSTRUMENT2
5150  double deta;
5151  if (_tiles[bj->tile_index].eta_centre == tile->eta_centre) deta = 0;
5152  else deta = std::abs(bj->eta - tile->eta_centre) - _tile_half_size_eta;
5153  double dphi = std::abs(bj->phi - tile->phi_centre);
5154  if (dphi > pi) dphi = twopi-dphi;
5155  dphi -= _tile_half_size_phi;
5156  if (dphi < 0) dphi = 0;
5157  return dphi*dphi + deta*deta;
5158 }
5159 inline void LazyTiling25::_update_jetX_jetI_NN(TiledJet * jetX, TiledJet * jetI, vector<TiledJet *> & jets_for_minheap) {
5160  double dist = _bj_dist(jetI,jetX);
5161  if (dist < jetI->NN_dist) {
5162  if (jetI != jetX) {
5163  jetI->NN_dist = dist;
5164  jetI->NN = jetX;
5165  if (!jetI->minheap_update_needed()) {
5167  jets_for_minheap.push_back(jetI);
5168  }
5169  }
5170  }
5171  if (dist < jetX->NN_dist) {
5172  if (jetI != jetX) {
5173  jetX->NN_dist = dist;
5174  jetX->NN = jetI;}
5175  }
5176 }
5177 inline void LazyTiling25::_set_NN(TiledJet * jetI,
5178  vector<TiledJet *> & jets_for_minheap) {
5179  jetI->NN_dist = _R2;
5180  jetI->NN = NULL;
5181  if (!jetI->minheap_update_needed()) {
5183  jets_for_minheap.push_back(jetI);}
5184  Tile25 * tile_ptr = &_tiles[jetI->tile_index];
5185  for (Tile25 ** near_tile = tile_ptr->begin_tiles;
5186  near_tile != tile_ptr->end_tiles; near_tile++) {
5187  if (jetI->NN_dist < _distance_to_tile(jetI, *near_tile)) continue;
5188  for (TiledJet * jetJ = (*near_tile)->head;
5189  jetJ != NULL; jetJ = jetJ->next) {
5190  double dist = _bj_dist(jetI,jetJ);
5191  if (dist < jetI->NN_dist && jetJ != jetI) {
5192  jetI->NN_dist = dist; jetI->NN = jetJ;
5193  }
5194  }
5195  }
5196 }
5198  int n = _jets.size();
5199  if (n == 0) return;
5200  TiledJet * briefjets = new TiledJet[n];
5201  TiledJet * jetA = briefjets, * jetB;
5202  TiledJet oldB = briefjets[0];
5203  vector<int> tile_union(3*25);
5204  for (int i = 0; i< n; i++) {
5205  _tj_set_jetinfo(jetA, i);
5206  jetA++; // move on to next entry of briefjets
5207  }
5208  TiledJet * head = briefjets; // a nicer way of naming start
5209  vector<Tile25>::iterator tile;
5210  for (tile = _tiles.begin(); tile != _tiles.end(); tile++) {
5211  for (jetA = tile->head; jetA != NULL; jetA = jetA->next) {
5212  for (jetB = tile->head; jetB != jetA; jetB = jetB->next) {
5213  double dist = _bj_dist_not_periodic(jetA,jetB);
5214  if (dist < jetA->NN_dist) {jetA->NN_dist = dist; jetA->NN = jetB;}
5215  if (dist < jetB->NN_dist) {jetB->NN_dist = dist; jetB->NN = jetA;}
5216  }
5217  }
5218  for (jetA = tile->head; jetA != NULL; jetA = jetA->next) {
5219  if (jetA->NN_dist > tile->max_NN_dist) tile->max_NN_dist = jetA->NN_dist;
5220  }
5221  }
5222  for (tile = _tiles.begin(); tile != _tiles.end(); tile++) {
5223  if (tile->use_periodic_delta_phi) {
5224  for (Tile25 ** RTile = tile->RH_tiles; RTile != tile->end_tiles; RTile++) {
5225  for (jetA = tile->head; jetA != NULL; jetA = jetA->next) {
5226  double dist_to_tile = _distance_to_tile(jetA, *RTile);
5227  bool relevant_for_jetA = dist_to_tile <= jetA->NN_dist;
5228  bool relevant_for_RTile = dist_to_tile <= (*RTile)->max_NN_dist;
5229  if (relevant_for_jetA || relevant_for_RTile) {
5230  for (jetB = (*RTile)->head; jetB != NULL; jetB = jetB->next) {
5231  double dist = _bj_dist(jetA,jetB);
5232  if (dist < jetA->NN_dist) {jetA->NN_dist = dist; jetA->NN = jetB;}
5233  if (dist < jetB->NN_dist) {jetB->NN_dist = dist; jetB->NN = jetA;}
5234  }
5235  }
5236  }
5237  }
5238  } else {
5239  for (Tile25 ** RTile = tile->RH_tiles; RTile != tile->end_tiles; RTile++) {
5240  for (jetA = tile->head; jetA != NULL; jetA = jetA->next) {
5241  double dist_to_tile = _distance_to_tile(jetA, *RTile);
5242  bool relevant_for_jetA = dist_to_tile <= jetA->NN_dist;
5243  bool relevant_for_RTile = dist_to_tile <= (*RTile)->max_NN_dist;
5244  if (relevant_for_jetA || relevant_for_RTile) {
5245  for (jetB = (*RTile)->head; jetB != NULL; jetB = jetB->next) {
5246  double dist = _bj_dist_not_periodic(jetA,jetB);
5247  if (dist < jetA->NN_dist) {jetA->NN_dist = dist; jetA->NN = jetB;}
5248  if (dist < jetB->NN_dist) {jetB->NN_dist = dist; jetB->NN = jetA;}
5249  }
5250  }
5251  }
5252  }
5253  }
5254  }
5255  for (tile = _tiles.begin(); tile != _tiles.end(); tile++) {
5256  tile->max_NN_dist = 0;
5257  for (jetA = tile->head; jetA != NULL; jetA = jetA->next) {
5258  if (jetA->NN_dist > tile->max_NN_dist) tile->max_NN_dist = jetA->NN_dist;
5259  }
5260  }
5261 #ifdef INSTRUMENT2
5262  cout << "intermediate ncall, dtt = " << _ncall << " " << _ncall_dtt << endl; // GPS tmp
5263 #endif // INSTRUMENT2
5264  vector<double> diJs(n);
5265  for (int i = 0; i < n; i++) {
5266  diJs[i] = _bj_diJ(&briefjets[i]);
5267  briefjets[i].label_minheap_update_done();
5268  }
5269  MinHeap minheap(diJs);
5270  vector<TiledJet *> jets_for_minheap;
5271  jets_for_minheap.reserve(n);
5272  int history_location = n-1;
5273  while (n > 0) {
5274  double diJ_min = minheap.minval() *_invR2;
5275  jetA = head + minheap.minloc();
5276  history_location++;
5277  jetB = jetA->NN;
5278  if (jetB != NULL) {
5279  if (jetA < jetB) {std::swap(jetA,jetB);}
5280  int nn; // new jet index
5281  _cs.plugin_record_ij_recombination(jetA->_jets_index, jetB->_jets_index, diJ_min, nn);
5282  _bj_remove_from_tiles(jetA);
5283  oldB = * jetB; // take a copy because we will need it...
5284  _bj_remove_from_tiles(jetB);
5285  _tj_set_jetinfo(jetB, nn); // cause jetB to become _jets[nn]
5286  } else {
5287  _cs.plugin_record_iB_recombination(jetA->_jets_index, diJ_min);
5288  _bj_remove_from_tiles(jetA);
5289  }
5290  minheap.remove(jetA-head);
5291  int n_near_tiles = 0;
5292  if (jetB != NULL) {
5293  Tile25 & jetB_tile = _tiles[jetB->tile_index];
5294  for (Tile25 ** near_tile = jetB_tile.begin_tiles;
5295  near_tile != jetB_tile.end_tiles; near_tile++) {
5296  double dist_to_tile = _distance_to_tile(jetB, *near_tile);
5297  bool relevant_for_jetB = dist_to_tile <= jetB->NN_dist;
5298  bool relevant_for_near_tile = dist_to_tile <= (*near_tile)->max_NN_dist;
5299  bool relevant = relevant_for_jetB || relevant_for_near_tile;
5300  if (! relevant) continue;
5301  tile_union[n_near_tiles] = *near_tile - & _tiles[0];
5302  (*near_tile)->tagged = true;
5303  n_near_tiles++;
5304  for (TiledJet * jetI = (*near_tile)->head; jetI != NULL; jetI = jetI->next) {
5305  if (jetI->NN == jetA || jetI->NN == jetB) _set_NN(jetI, jets_for_minheap);
5306  _update_jetX_jetI_NN(jetB, jetI, jets_for_minheap);
5307  }
5308  }
5309  }
5310  int n_done_tiles = n_near_tiles;
5312  tile_union, n_near_tiles);
5313  if (jetB != NULL) {
5315  tile_union,n_near_tiles);
5316  jetB->label_minheap_update_needed();
5317  jets_for_minheap.push_back(jetB);
5318  }
5319  for (int itile = 0; itile < n_done_tiles; itile++) {
5320  _tiles[tile_union[itile]].tagged = false;
5321  }
5322  for (int itile = n_done_tiles; itile < n_near_tiles; itile++) {
5323  Tile25 * tile_ptr = &_tiles[tile_union[itile]];
5324  tile_ptr->tagged = false;
5325  for (TiledJet * jetI = tile_ptr->head; jetI != NULL; jetI = jetI->next) {
5326  if (jetI->NN == jetA || (jetI->NN == jetB && jetB != NULL)) {
5327  _set_NN(jetI, jets_for_minheap);
5328  }
5329  }
5330  }
5331  while (jets_for_minheap.size() > 0) {
5332  TiledJet * jetI = jets_for_minheap.back();
5333  jets_for_minheap.pop_back();
5334  minheap.update(jetI-head, _bj_diJ(jetI));
5335  jetI->label_minheap_update_done();
5336  Tile25 & tile_I = _tiles[jetI->tile_index];
5337  if (tile_I.max_NN_dist < jetI->NN_dist) tile_I.max_NN_dist = jetI->NN_dist;
5338  }
5339  n--;
5340  }
5341  delete[] briefjets;
5342 #ifdef INSTRUMENT2
5343  cout << "ncall, dtt = " << _ncall << " " << _ncall_dtt << endl; // GPS tmp
5344 #endif // INSTRUMENT2
5345 }
5347 #include <iomanip>
5348 #include <limits>
5349 #include <cmath>
5350 using namespace std;
5351 #define _FJCORE_TILING2_USE_TILING_ANALYSIS_
5352 FJCORE_BEGIN_NAMESPACE // defined in fastjet/internal/base.hh
5353 LazyTiling9::LazyTiling9(ClusterSequence & cs) :
5354  _cs(cs), _jets(cs.jets())
5355 {
5356 #ifdef INSTRUMENT2
5357  _ncall = 0; // gps tmp
5358  _ncall_dtt = 0; // gps tmp
5359 #endif // INSTRUMENT2
5360  _Rparam = cs.jet_def().R();
5361  _R2 = _Rparam * _Rparam;
5362  _invR2 = 1.0 / _R2;
5364 }
5366  double default_size = max(0.1,_Rparam);
5367  _tile_size_eta = default_size;
5368  _n_tiles_phi = max(3,int(floor(twopi/default_size)));
5369  _tile_size_phi = twopi / _n_tiles_phi; // >= _Rparam and fits in 2pi
5370 #ifdef _FJCORE_TILING2_USE_TILING_ANALYSIS_
5371  TilingExtent tiling_analysis(_cs);
5372  _tiles_eta_min = tiling_analysis.minrap();
5373  _tiles_eta_max = tiling_analysis.maxrap();
5374 #else
5375  _tiles_eta_min = 0.0;
5376  _tiles_eta_max = 0.0;
5377  const double maxrap = 7.0;
5378  for(unsigned int i = 0; i < _jets.size(); i++) {
5379  double eta = _jets[i].rap();
5380  if (abs(eta) < maxrap) {
5381  if (eta < _tiles_eta_min) {_tiles_eta_min = eta;}
5382  if (eta > _tiles_eta_max) {_tiles_eta_max = eta;}
5383  }
5384  }
5385 #endif
5386 # define FJCORE_LAZY9_MIN2TILESY
5387 #ifdef FJCORE_LAZY9_MIN2TILESY
5390  _tiles_ieta_min = 0;
5391  _tiles_ieta_max = 1;
5393  } else {
5394 #endif //FASTJET_LAZY9_MIN2TILESY
5399 #ifdef FJCORE_LAZY9_MIN2TILESY
5400  }
5401 #endif
5404  vector<bool> use_periodic_delta_phi(_n_tiles_phi, false);
5405  if (_n_tiles_phi <= 3) {
5406  fill(use_periodic_delta_phi.begin(), use_periodic_delta_phi.end(), true);
5407  } else {
5408  use_periodic_delta_phi[0] = true;
5409  use_periodic_delta_phi[_n_tiles_phi-1] = true;
5410  }
5412  for (int ieta = _tiles_ieta_min; ieta <= _tiles_ieta_max; ieta++) {
5413  for (int iphi = 0; iphi < _n_tiles_phi; iphi++) {
5414  Tile2 * tile = & _tiles[_tile_index(ieta,iphi)];
5415  tile->head = NULL; // first element of tiles points to itself
5416  tile->begin_tiles[0] = tile;
5417  Tile2 ** pptile = & (tile->begin_tiles[0]);
5418  pptile++;
5419  tile->surrounding_tiles = pptile;
5420  if (ieta > _tiles_ieta_min) {
5421  // with the itile subroutine, we can safely run tiles from
5422  // idphi=-1 to idphi=+1, because it takes care of
5423  // negative and positive boundaries
5424  for (int idphi = -1; idphi <=+1; idphi++) {
5425  *pptile = & _tiles[_tile_index(ieta-1,iphi+idphi)];
5426  pptile++;
5427  }
5428  }
5429  *pptile = & _tiles[_tile_index(ieta,iphi-1)];
5430  pptile++;
5431  tile->RH_tiles = pptile;
5432  *pptile = & _tiles[_tile_index(ieta,iphi+1)];
5433  pptile++;
5434  if (ieta < _tiles_ieta_max) {
5435  for (int idphi = -1; idphi <= +1; idphi++) {
5436  *pptile = & _tiles[_tile_index(ieta+1,iphi+idphi)];
5437  pptile++;
5438  }
5439  }
5440  tile->end_tiles = pptile;
5441  tile->tagged = false;
5442  tile->use_periodic_delta_phi = use_periodic_delta_phi[iphi];
5443  tile->max_NN_dist = 0;
5445  tile->phi_centre = (iphi+0.5)*_tile_size_phi;
5446  }
5447  }
5448 }
5449 int LazyTiling9::_tile_index(const double eta, const double phi) const {
5450  int ieta, iphi;
5451  if (eta <= _tiles_eta_min) {ieta = 0;}
5452  else if (eta >= _tiles_eta_max) {ieta = _tiles_ieta_max-_tiles_ieta_min;}
5453  else {
5454  ieta = int(((eta - _tiles_eta_min) / _tile_size_eta));
5455  if (ieta > _tiles_ieta_max-_tiles_ieta_min) {
5457  }
5458  iphi = int((phi+twopi)/_tile_size_phi) % _n_tiles_phi;
5459  return (iphi + ieta * _n_tiles_phi);
5460 }
5461 inline void LazyTiling9::_tj_set_jetinfo( TiledJet * const jet,
5462  const int _jets_index) {
5463  _bj_set_jetinfo<>(jet, _jets_index);
5464  jet->tile_index = _tile_index(jet->eta, jet->phi);
5465  Tile2 * tile = &_tiles[jet->tile_index];
5466  jet->previous = NULL;
5467  jet->next = tile->head;
5468  if (jet->next != NULL) {jet->next->previous = jet;}
5469  tile->head = jet;
5470 }
5472  Tile2 * tile = & _tiles[jet->tile_index];
5473  if (jet->previous == NULL) {
5474  tile->head = jet->next;
5475  } else {
5476  jet->previous->next = jet->next;
5477  }
5478  if (jet->next != NULL) {
5479  jet->next->previous = jet->previous;
5480  }
5481 }
5482 void LazyTiling9::_print_tiles(TiledJet * briefjets ) const {
5483  for (vector<Tile2>::const_iterator tile = _tiles.begin();
5484  tile < _tiles.end(); tile++) {
5485  cout << "Tile " << tile - _tiles.begin()<<" = ";
5486  vector<int> list;
5487  for (TiledJet * jetI = tile->head; jetI != NULL; jetI = jetI->next) {
5488  list.push_back(jetI-briefjets);
5489  }
5490  sort(list.begin(),list.end());
5491  for (unsigned int i = 0; i < list.size(); i++) {cout <<" "<<list[i];}
5492  cout <<"\n";
5493  }
5494 }
5496  vector<int> & tile_union, int & n_near_tiles) const {
5497  for (Tile2 * const * near_tile = _tiles[tile_index].begin_tiles;
5498  near_tile != _tiles[tile_index].end_tiles; near_tile++){
5499  tile_union[n_near_tiles] = *near_tile - & _tiles[0];
5500  n_near_tiles++;
5501  }
5502 }
5504  const int tile_index,
5505  vector<int> & tile_union, int & n_near_tiles) {
5506  for (Tile2 ** near_tile = _tiles[tile_index].begin_tiles;
5507  near_tile != _tiles[tile_index].end_tiles; near_tile++){
5508  if (! (*near_tile)->tagged) {
5509  (*near_tile)->tagged = true;
5510  tile_union[n_near_tiles] = *near_tile - & _tiles[0];
5511  n_near_tiles++;
5512  }
5513  }
5514 }
5516  const TiledJet * jet,
5517  vector<int> & tile_union, int & n_near_tiles) {
5518  Tile2 & tile = _tiles[jet->tile_index];
5519  for (Tile2 ** near_tile = tile.begin_tiles; near_tile != tile.end_tiles; near_tile++){
5520  if ((*near_tile)->tagged) continue;
5521  double dist = _distance_to_tile(jet, *near_tile) - tile_edge_security_margin;
5522  if (dist > (*near_tile)->max_NN_dist) continue;
5523  (*near_tile)->tagged = true;
5524  tile_union[n_near_tiles] = *near_tile - & _tiles[0];
5525  n_near_tiles++;
5526  }
5527 }
5528 inline double LazyTiling9::_distance_to_tile(const TiledJet * bj, const Tile2 * tile)
5529 #ifdef INSTRUMENT2
5530  {
5531  _ncall_dtt++; // GPS tmp
5532 #else
5533  const {
5534 #endif // INSTRUMENT2
5535  double deta;
5536  if (_tiles[bj->tile_index].eta_centre == tile->eta_centre) deta = 0;
5537  else deta = std::abs(bj->eta - tile->eta_centre) - _tile_half_size_eta;
5538  double dphi = std::abs(bj->phi - tile->phi_centre);
5539  if (dphi > pi) dphi = twopi-dphi;
5540  dphi -= _tile_half_size_phi;
5541  if (dphi < 0) dphi = 0;
5542  return dphi*dphi + deta*deta;
5543 }
5544 inline void LazyTiling9::_update_jetX_jetI_NN(TiledJet * jetX, TiledJet * jetI, vector<TiledJet *> & jets_for_minheap) {
5545  double dist = _bj_dist(jetI,jetX);
5546  if (dist < jetI->NN_dist) {
5547  if (jetI != jetX) {
5548  jetI->NN_dist = dist;
5549  jetI->NN = jetX;
5550  if (!jetI->minheap_update_needed()) {
5552  jets_for_minheap.push_back(jetI);
5553  }
5554  }
5555  }
5556  if (dist < jetX->NN_dist) {
5557  if (jetI != jetX) {
5558  jetX->NN_dist = dist;
5559  jetX->NN = jetI;}
5560  }
5561 }
5562 inline void LazyTiling9::_set_NN(TiledJet * jetI,
5563  vector<TiledJet *> & jets_for_minheap) {
5564  jetI->NN_dist = _R2;
5565  jetI->NN = NULL;
5566  if (!jetI->minheap_update_needed()) {
5568  jets_for_minheap.push_back(jetI);}
5569  Tile2 * tile_ptr = &_tiles[jetI->tile_index];
5570  for (Tile2 ** near_tile = tile_ptr->begin_tiles;
5571  near_tile != tile_ptr->end_tiles; near_tile++) {
5572  if (jetI->NN_dist < _distance_to_tile(jetI, *near_tile)) continue;
5573  for (TiledJet * jetJ = (*near_tile)->head;
5574  jetJ != NULL; jetJ = jetJ->next) {
5575  double dist = _bj_dist(jetI,jetJ);
5576  if (dist < jetI->NN_dist && jetJ != jetI) {
5577  jetI->NN_dist = dist; jetI->NN = jetJ;
5578  }
5579  }
5580  }
5581 }
5583  int n = _jets.size();
5584  if (n == 0) return;
5585  TiledJet * briefjets = new TiledJet[n];
5586  TiledJet * jetA = briefjets, * jetB;
5587  TiledJet oldB = briefjets[0];
5588  vector<int> tile_union(3*n_tile_neighbours);
5589  for (int i = 0; i< n; i++) {
5590  _tj_set_jetinfo(jetA, i);
5591  jetA++; // move on to next entry of briefjets
5592  }
5593  TiledJet * head = briefjets; // a nicer way of naming start
5594  vector<Tile2>::iterator tile;
5595  for (tile = _tiles.begin(); tile != _tiles.end(); tile++) {
5596  for (jetA = tile->head; jetA != NULL; jetA = jetA->next) {
5597  for (jetB = tile->head; jetB != jetA; jetB = jetB->next) {
5598  double dist = _bj_dist_not_periodic(jetA,jetB);
5599  if (dist < jetA->NN_dist) {jetA->NN_dist = dist; jetA->NN = jetB;}
5600  if (dist < jetB->NN_dist) {jetB->NN_dist = dist; jetB->NN = jetA;}
5601  }
5602  }
5603  for (jetA = tile->head; jetA != NULL; jetA = jetA->next) {
5604  if (jetA->NN_dist > tile->max_NN_dist) tile->max_NN_dist = jetA->NN_dist;
5605  }
5606  }
5607  for (tile = _tiles.begin(); tile != _tiles.end(); tile++) {
5608  if (tile->use_periodic_delta_phi) {
5609  for (Tile2 ** RTile = tile->RH_tiles; RTile != tile->end_tiles; RTile++) {
5610  for (jetA = tile->head; jetA != NULL; jetA = jetA->next) {
5611  double dist_to_tile = _distance_to_tile(jetA, *RTile);
5612  bool relevant_for_jetA = dist_to_tile <= jetA->NN_dist;
5613  bool relevant_for_RTile = dist_to_tile <= (*RTile)->max_NN_dist;
5614  if (relevant_for_jetA || relevant_for_RTile) {
5615  for (jetB = (*RTile)->head; jetB != NULL; jetB = jetB->next) {
5616  double dist = _bj_dist(jetA,jetB);
5617  if (dist < jetA->NN_dist) {jetA->NN_dist = dist; jetA->NN = jetB;}
5618  if (dist < jetB->NN_dist) {jetB->NN_dist = dist; jetB->NN = jetA;}
5619  }
5620  }
5621  }
5622  }
5623  } else {
5624  for (Tile2 ** RTile = tile->RH_tiles; RTile != tile->end_tiles; RTile++) {
5625  for (jetA = tile->head; jetA != NULL; jetA = jetA->next) {
5626  double dist_to_tile = _distance_to_tile(jetA, *RTile);
5627  bool relevant_for_jetA = dist_to_tile <= jetA->NN_dist;
5628  bool relevant_for_RTile = dist_to_tile <= (*RTile)->max_NN_dist;
5629  if (relevant_for_jetA || relevant_for_RTile) {
5630  for (jetB = (*RTile)->head; jetB != NULL; jetB = jetB->next) {
5631  double dist = _bj_dist_not_periodic(jetA,jetB);
5632  if (dist < jetA->NN_dist) {jetA->NN_dist = dist; jetA->NN = jetB;}
5633  if (dist < jetB->NN_dist) {jetB->NN_dist = dist; jetB->NN = jetA;}
5634  }
5635  }
5636  }
5637  }
5638  }
5639  }
5640  for (tile = _tiles.begin(); tile != _tiles.end(); tile++) {
5641  tile->max_NN_dist = 0;
5642  for (jetA = tile->head; jetA != NULL; jetA = jetA->next) {
5643  if (jetA->NN_dist > tile->max_NN_dist) tile->max_NN_dist = jetA->NN_dist;
5644  }
5645  }
5646 #ifdef INSTRUMENT2
5647  cout << "intermediate ncall, dtt = " << _ncall << " " << _ncall_dtt << endl; // GPS tmp
5648 #endif // INSTRUMENT2
5649  vector<double> diJs(n);
5650  for (int i = 0; i < n; i++) {
5651  diJs[i] = _bj_diJ(&briefjets[i]);
5652  briefjets[i].label_minheap_update_done();
5653  }
5654  MinHeap minheap(diJs);
5655  vector<TiledJet *> jets_for_minheap;
5656  jets_for_minheap.reserve(n);
5657  int history_location = n-1;
5658  while (n > 0) {
5659  double diJ_min = minheap.minval() *_invR2;
5660  jetA = head + minheap.minloc();
5661  history_location++;
5662  jetB = jetA->NN;
5663  if (jetB != NULL) {
5664  if (jetA < jetB) {std::swap(jetA,jetB);}
5665  int nn; // new jet index
5666  _cs.plugin_record_ij_recombination(jetA->_jets_index, jetB->_jets_index, diJ_min, nn);
5667  _bj_remove_from_tiles(jetA);
5668  oldB = * jetB; // take a copy because we will need it...
5669  _bj_remove_from_tiles(jetB);
5670  _tj_set_jetinfo(jetB, nn); // cause jetB to become _jets[nn]
5671  } else {
5672  _cs.plugin_record_iB_recombination(jetA->_jets_index, diJ_min);
5673  _bj_remove_from_tiles(jetA);
5674  }
5675  minheap.remove(jetA-head);
5676  int n_near_tiles = 0;
5677  if (jetB != NULL) {
5678  Tile2 & jetB_tile = _tiles[jetB->tile_index];
5679  for (Tile2 ** near_tile = jetB_tile.begin_tiles;
5680  near_tile != jetB_tile.end_tiles; near_tile++) {
5681  double dist_to_tile = _distance_to_tile(jetB, *near_tile);
5682  bool relevant_for_jetB = dist_to_tile <= jetB->NN_dist;
5683  bool relevant_for_near_tile = dist_to_tile <= (*near_tile)->max_NN_dist;
5684  bool relevant = relevant_for_jetB || relevant_for_near_tile;
5685  if (! relevant) continue;
5686  tile_union[n_near_tiles] = *near_tile - & _tiles[0];
5687  (*near_tile)->tagged = true;
5688  n_near_tiles++;
5689  for (TiledJet * jetI = (*near_tile)->head; jetI != NULL; jetI = jetI->next) {
5690  if (jetI->NN == jetA || jetI->NN == jetB) _set_NN(jetI, jets_for_minheap);
5691  _update_jetX_jetI_NN(jetB, jetI, jets_for_minheap);
5692  }
5693  }
5694  }
5695  int n_done_tiles = n_near_tiles;
5697  tile_union, n_near_tiles);
5698  if (jetB != NULL) {
5700  tile_union,n_near_tiles);
5701  jetB->label_minheap_update_needed();
5702  jets_for_minheap.push_back(jetB);
5703  }
5704  for (int itile = 0; itile < n_done_tiles; itile++) {
5705  _tiles[tile_union[itile]].tagged = false;
5706  }
5707  for (int itile = n_done_tiles; itile < n_near_tiles; itile++) {
5708  Tile2 * tile_ptr = &_tiles[tile_union[itile]];
5709  tile_ptr->tagged = false;
5710  for (TiledJet * jetI = tile_ptr->head; jetI != NULL; jetI = jetI->next) {
5711  if (jetI->NN == jetA || (jetI->NN == jetB && jetB != NULL)) {
5712  _set_NN(jetI, jets_for_minheap);
5713  }
5714  }
5715  }
5716  while (jets_for_minheap.size() > 0) {
5717  TiledJet * jetI = jets_for_minheap.back();
5718  jets_for_minheap.pop_back();
5719  minheap.update(jetI-head, _bj_diJ(jetI));
5720  jetI->label_minheap_update_done();
5721  Tile2 & tile_I = _tiles[jetI->tile_index];
5722  if (tile_I.max_NN_dist < jetI->NN_dist) tile_I.max_NN_dist = jetI->NN_dist;
5723  }
5724  n--;
5725  }
5726  delete[] briefjets;
5727 #ifdef INSTRUMENT2
5728  cout << "ncall, dtt = " << _ncall << " " << _ncall_dtt << endl; // GPS tmp
5729 #endif // INSTRUMENT2
5730 }
5732 #include <iomanip>
5733 using namespace std;
5734 FJCORE_BEGIN_NAMESPACE // defined in fastjet/internal/base.hh
5735 LazyTiling9Alt::LazyTiling9Alt(ClusterSequence & cs) :
5736  _cs(cs), _jets(cs.jets())
5737 {
5738  _Rparam = cs.jet_def().R();
5739  _R2 = _Rparam * _Rparam;
5740  _invR2 = 1.0 / _R2;
5742 }
5744  double default_size = max(0.1,_Rparam);
5745  _tile_size_eta = default_size;
5746  _n_tiles_phi = max(3,int(floor(twopi/default_size)));
5747  _tile_size_phi = twopi / _n_tiles_phi; // >= _Rparam and fits in 2pi
5748  _tiles_eta_min = 0.0;
5749  _tiles_eta_max = 0.0;
5750  const double maxrap = 7.0;
5751  for(unsigned int i = 0; i < _jets.size(); i++) {
5752  double eta = _jets[i].rap();
5753  if (abs(eta) < maxrap) {
5754  if (eta < _tiles_eta_min) {_tiles_eta_min = eta;}
5755  if (eta > _tiles_eta_max) {_tiles_eta_max = eta;}
5756  }
5757  }
5762  _tile_half_size_eta = _tile_size_eta * 0.5;
5764  vector<bool> use_periodic_delta_phi(_n_tiles_phi, false);
5765  if (_n_tiles_phi <= 3) {
5766  fill(use_periodic_delta_phi.begin(), use_periodic_delta_phi.end(), true);
5767  } else {
5768  use_periodic_delta_phi[0] = true;
5769  use_periodic_delta_phi[_n_tiles_phi-1] = true;
5770  }
5772  for (int ieta = _tiles_ieta_min; ieta <= _tiles_ieta_max; ieta++) {
5773  for (int iphi = 0; iphi < _n_tiles_phi; iphi++) {
5774  Tile * tile = & _tiles[_tile_index(ieta,iphi)];
5775  tile->head = NULL; // first element of tiles points to itself
5777  Tile::TileFnPair * pptile = & (tile->begin_tiles[0]);
5778  pptile++;
5779  tile->surrounding_tiles = pptile;
5780  if (ieta > _tiles_ieta_min) {
5781  // with the itile subroutine, we can safely run tiles from
5782  // idphi=-1 to idphi=+1, because it takes care of
5783  // negative and positive boundaries
5784  //for (int idphi = -1; idphi <=+1; idphi++) {
5785  *pptile = Tile::TileFnPair(& _tiles[_tile_index(ieta-1,iphi-1)],
5787  pptile++;
5788  *pptile = Tile::TileFnPair(& _tiles[_tile_index(ieta-1,iphi)],
5790  pptile++;
5791  *pptile = Tile::TileFnPair(& _tiles[_tile_index(ieta-1,iphi+1)],
5793  pptile++;
5794  }
5795  *pptile = Tile::TileFnPair(& _tiles[_tile_index(ieta,iphi-1)],
5797  pptile++;
5798  tile->RH_tiles = pptile;
5799  *pptile = Tile::TileFnPair(& _tiles[_tile_index(ieta,iphi+1)],
5801  pptile++;
5802  if (ieta < _tiles_ieta_max) {
5803  //for (int idphi = -1; idphi <= +1; idphi++) {
5804  // *pptile = & _tiles[_tile_index(ieta+1,iphi+idphi)];
5805  // pptile++;
5806  //}
5807  *pptile = Tile::TileFnPair(& _tiles[_tile_index(ieta+1,iphi-1)],
5809  pptile++;
5810  *pptile = Tile::TileFnPair(& _tiles[_tile_index(ieta+1,iphi)],
5812  pptile++;
5813  *pptile = Tile::TileFnPair(& _tiles[_tile_index(ieta+1,iphi+1)],
5815  pptile++;
5816  }
5817  tile->end_tiles = pptile;
5818  tile->tagged = false;
5819  tile->use_periodic_delta_phi = use_periodic_delta_phi[iphi];
5820  tile->max_NN_dist = 0;
5821  tile->eta_min = ieta*_tile_size_eta;
5822  tile->eta_max = (ieta+1)*_tile_size_eta;
5823  tile->phi_min = iphi*_tile_size_phi;
5824  tile->phi_max = (iphi+1)*_tile_size_phi;
5825  }
5826  }
5827 }
5828 int LazyTiling9Alt::_tile_index(const double eta, const double phi) const {
5829  int ieta, iphi;
5830  if (eta <= _tiles_eta_min) {ieta = 0;}
5831  else if (eta >= _tiles_eta_max) {ieta = _tiles_ieta_max-_tiles_ieta_min;}
5832  else {
5833  ieta = int(((eta - _tiles_eta_min) / _tile_size_eta));
5834  if (ieta > _tiles_ieta_max-_tiles_ieta_min) {
5836  }
5837  iphi = int((phi+twopi)/_tile_size_phi) % _n_tiles_phi;
5838  return (iphi + ieta * _n_tiles_phi);
5839 }
5841  const int _jets_index) {
5842  _bj_set_jetinfo<>(jet, _jets_index);
5843  jet->tile_index = _tile_index(jet->eta, jet->phi);
5844  Tile * tile = &_tiles[jet->tile_index];
5845  jet->previous = NULL;
5846  jet->next = tile->head;
5847  if (jet->next != NULL) {jet->next->previous = jet;}
5848  tile->head = jet;
5849 }
5851  Tile * tile = & _tiles[jet->tile_index];
5852  if (jet->previous == NULL) {
5853  tile->head = jet->next;
5854  } else {
5855  jet->previous->next = jet->next;
5856  }
5857  if (jet->next != NULL) {
5858  jet->next->previous = jet->previous;
5859  }
5860 }
5861 void LazyTiling9Alt::_print_tiles(TiledJet * briefjets ) const {
5862  for (vector<Tile>::const_iterator tile = _tiles.begin();
5863  tile < _tiles.end(); tile++) {
5864  cout << "Tile " << tile - _tiles.begin()<<" = ";
5865  vector<int> list;
5866  for (TiledJet * jetI = tile->head; jetI != NULL; jetI = jetI->next) {
5867  list.push_back(jetI-briefjets);
5868  }
5869  sort(list.begin(),list.end());
5870  for (unsigned int i = 0; i < list.size(); i++) {cout <<" "<<list[i];}
5871  cout <<"\n";
5872  }
5873 }
5875  vector<int> & tile_union, int & n_near_tiles) const {
5876  for (Tile::TileFnPair const * near_tile = _tiles[tile_index].begin_tiles;
5877  near_tile != _tiles[tile_index].end_tiles; near_tile++){
5878  tile_union[n_near_tiles] = near_tile->first - & _tiles[0];
5879  n_near_tiles++;
5880  }
5881 }
5883  const int tile_index,
5884  vector<int> & tile_union, int & n_near_tiles) {
5885  for (Tile::TileFnPair * near_tile = _tiles[tile_index].begin_tiles;
5886  near_tile != _tiles[tile_index].end_tiles; near_tile++){
5887  if (! (near_tile->first)->tagged) {
5888  (near_tile->first)->tagged = true;
5889  tile_union[n_near_tiles] = near_tile->first - & _tiles[0];
5890  n_near_tiles++;
5891  }
5892  }
5893 }
5895  const TiledJet * jet,
5896  vector<int> & tile_union, int & n_near_tiles) {
5897  Tile & tile = _tiles[jet->tile_index];
5898  for (Tile::TileFnPair * near_tile = tile.begin_tiles; near_tile != tile.end_tiles; near_tile++){
5899  if ((near_tile->first)->tagged) continue;
5900  double dist = (tile.*(near_tile->second))(jet) - tile_edge_security_margin;
5901  if (dist > (near_tile->first)->max_NN_dist) continue;
5902  (near_tile->first)->tagged = true;
5903  tile_union[n_near_tiles] = near_tile->first - & _tiles[0];
5904  n_near_tiles++;
5905  }
5906 }
5907 ostream & operator<<(ostream & ostr, const TiledJet & jet) {
5908  ostr << "j" << setw(3) << jet._jets_index << ":pt2,rap,phi=" ; ostr.flush();
5909  ostr << jet.kt2 << ","; ostr.flush();
5910  ostr << jet.eta << ","; ostr.flush();
5911  ostr << jet.phi; ostr.flush();
5912  ostr << ", tile=" << jet.tile_index; ostr.flush();
5913  return ostr;
5914 }
5915 inline void LazyTiling9Alt::_update_jetX_jetI_NN(TiledJet * jetX, TiledJet * jetI, vector<TiledJet *> & jets_for_minheap) {
5916  double dist = _bj_dist(jetI,jetX);
5917  if (dist < jetI->NN_dist) {
5918  if (jetI != jetX) {
5919  jetI->NN_dist = dist;
5920  jetI->NN = jetX;
5921  if (!jetI->minheap_update_needed()) {
5923  jets_for_minheap.push_back(jetI);
5924  }
5925  }
5926  }
5927  if (dist < jetX->NN_dist) {
5928  if (jetI != jetX) {
5929  jetX->NN_dist = dist;
5930  jetX->NN = jetI;}
5931  }
5932 }
5933 inline void LazyTiling9Alt::_set_NN(TiledJet * jetI,
5934  vector<TiledJet *> & jets_for_minheap) {
5935  jetI->NN_dist = _R2;
5936  jetI->NN = NULL;
5937  if (!jetI->minheap_update_needed()) {
5939  jets_for_minheap.push_back(jetI);}
5940  Tile * tile_ptr = &_tiles[jetI->tile_index];
5941  for (Tile::TileFnPair * near_tile = tile_ptr->begin_tiles;
5942  near_tile != tile_ptr->end_tiles; near_tile++) {
5943  if (jetI->NN_dist < (tile_ptr->*(near_tile->second))(jetI)) continue;
5944  for (TiledJet * jetJ = (near_tile->first)->head;
5945  jetJ != NULL; jetJ = jetJ->next) {
5946  double dist = _bj_dist(jetI,jetJ);
5947  if (dist < jetI->NN_dist && jetJ != jetI) {
5948  jetI->NN_dist = dist; jetI->NN = jetJ;
5949  }
5950  }
5951  }
5952 }
5954  int n = _jets.size();
5955  TiledJet * briefjets = new TiledJet[n];
5956  TiledJet * jetA = briefjets, * jetB;
5957  TiledJet oldB;
5958  vector<int> tile_union(3*n_tile_neighbours);
5959  for (int i = 0; i< n; i++) {
5960  _tj_set_jetinfo(jetA, i);
5961  jetA++; // move on to next entry of briefjets
5962  }
5963  TiledJet * head = briefjets; // a nicer way of naming start
5964  vector<Tile>::iterator tile;
5965  for (tile = _tiles.begin(); tile != _tiles.end(); tile++) {
5966  for (jetA = tile->head; jetA != NULL; jetA = jetA->next) {
5967  for (jetB = tile->head; jetB != jetA; jetB = jetB->next) {
5968  double dist = _bj_dist_not_periodic(jetA,jetB);
5969  if (dist < jetA->NN_dist) {jetA->NN_dist = dist; jetA->NN = jetB;}
5970  if (dist < jetB->NN_dist) {jetB->NN_dist = dist; jetB->NN = jetA;}
5971  }
5972  }
5973  for (jetA = tile->head; jetA != NULL; jetA = jetA->next) {
5974  if (jetA->NN_dist > tile->max_NN_dist) tile->max_NN_dist = jetA->NN_dist;
5975  }
5976  }
5977  for (tile = _tiles.begin(); tile != _tiles.end(); tile++) {
5978  if (tile->use_periodic_delta_phi) {
5979  for (Tile::TileFnPair * RTileFnPair = tile->RH_tiles;
5980  RTileFnPair != tile->end_tiles; RTileFnPair++) {
5981  Tile *RTile = RTileFnPair->first;
5982  for (jetA = tile->head; jetA != NULL; jetA = jetA->next) {
5983  double dist_to_tile = ((*tile).*(RTileFnPair->second))(jetA);
5984  bool relevant_for_jetA = dist_to_tile <= jetA->NN_dist;
5985  bool relevant_for_RTile = dist_to_tile <= RTile->max_NN_dist;
5986  if (relevant_for_jetA || relevant_for_RTile) {
5987  for (jetB = RTile->head; jetB != NULL; jetB = jetB->next) {
5988  double dist = _bj_dist(jetA,jetB);
5989  if (dist < jetA->NN_dist) {jetA->NN_dist = dist; jetA->NN = jetB;}
5990  if (dist < jetB->NN_dist) {jetB->NN_dist = dist; jetB->NN = jetA;}
5991  }
5992  }
5993  }
5994  }
5995  } else {
5996  for (Tile::TileFnPair* RTileFnPair = tile->RH_tiles;
5997  RTileFnPair != tile->end_tiles; RTileFnPair++) {
5998  Tile *RTile = RTileFnPair->first;
5999  for (jetA = tile->head; jetA != NULL; jetA = jetA->next) {
6000  double dist_to_tile = ((*tile).*(RTileFnPair->second))(jetA);
6001  bool relevant_for_jetA = dist_to_tile <= jetA->NN_dist;
6002  bool relevant_for_RTile = dist_to_tile <= RTile->max_NN_dist;
6003  if (relevant_for_jetA || relevant_for_RTile) {
6004  for (jetB = RTile->head; jetB != NULL; jetB = jetB->next) {
6005  double dist = _bj_dist_not_periodic(jetA,jetB);
6006  if (dist < jetA->NN_dist) {jetA->NN_dist = dist; jetA->NN = jetB;}
6007  if (dist < jetB->NN_dist) {jetB->NN_dist = dist; jetB->NN = jetA;}
6008  }
6009  }
6010  }
6011  }
6012  }
6013  }
6014  for (tile = _tiles.begin(); tile != _tiles.end(); tile++) {
6015  tile->max_NN_dist = 0;
6016  for (jetA = tile->head; jetA != NULL; jetA = jetA->next) {
6017  if (jetA->NN_dist > tile->max_NN_dist) tile->max_NN_dist = jetA->NN_dist;
6018  }
6019  }
6020  vector<double> diJs(n);
6021  for (int i = 0; i < n; i++) {
6022  diJs[i] = _bj_diJ(&briefjets[i]);
6023  briefjets[i].label_minheap_update_done();
6024  }
6025  MinHeap minheap(diJs);
6026  vector<TiledJet *> jets_for_minheap;
6027  jets_for_minheap.reserve(n);
6028  int history_location = n-1;
6029  while (n > 0) {
6030  double diJ_min = minheap.minval() *_invR2;
6031  jetA = head + minheap.minloc();
6032  history_location++;
6033  jetB = jetA->NN;
6034  if (jetB != NULL) {
6035  if (jetA < jetB) {std::swap(jetA,jetB);}
6036  int nn; // new jet index
6037  _cs.plugin_record_ij_recombination(jetA->_jets_index, jetB->_jets_index, diJ_min, nn);
6038  _bj_remove_from_tiles(jetA);
6039  oldB = * jetB; // take a copy because we will need it...
6040  _bj_remove_from_tiles(jetB);
6041  _tj_set_jetinfo(jetB, nn); // cause jetB to become _jets[nn]
6042  } else {
6043  _cs.plugin_record_iB_recombination(jetA->_jets_index, diJ_min);
6044  _bj_remove_from_tiles(jetA);
6045  }
6046  minheap.remove(jetA-head);
6047  int n_near_tiles = 0;
6049  tile_union, n_near_tiles);
6050  if (jetB != NULL) {
6052  tile_union,n_near_tiles);
6053  jetB->label_minheap_update_needed();
6054  jets_for_minheap.push_back(jetB);
6055  }
6056  if (jetB != NULL) {
6057  Tile & jetB_tile = _tiles[jetB->tile_index];
6058  for (Tile::TileFnPair * near_tile_fn_pair = jetB_tile.begin_tiles;
6059  near_tile_fn_pair != jetB_tile.end_tiles; near_tile_fn_pair++) {
6060  Tile * near_tile = near_tile_fn_pair->first;
6061  double dist_to_tile = (jetB_tile.*(near_tile_fn_pair->second))(jetB);
6062  bool relevant_for_jetB = dist_to_tile <= jetB->NN_dist;
6063  bool relevant_for_near_tile = dist_to_tile <= near_tile->max_NN_dist;
6064  bool relevant = relevant_for_jetB || relevant_for_near_tile;
6065  if (relevant) {
6066  if (near_tile->tagged) {
6067  for (TiledJet * jetI = near_tile->head; jetI != NULL; jetI = jetI->next) {
6068  if (jetI->NN == jetA || jetI->NN == jetB) _set_NN(jetI, jets_for_minheap);
6069  _update_jetX_jetI_NN(jetB, jetI, jets_for_minheap);
6070  }
6071  near_tile->tagged = false;
6072  } else {
6073  for (TiledJet * jetI = near_tile->head; jetI != NULL; jetI = jetI->next) {
6074  _update_jetX_jetI_NN(jetB, jetI, jets_for_minheap);
6075  }
6076  }
6077  }
6078  // // -- Keep this old inline code for later speed tests
6079  }
6080  }
6081  for (int itile = 0; itile < n_near_tiles; itile++) {
6082  Tile * tile_ptr = &_tiles[tile_union[itile]];
6083  if (!tile_ptr->tagged) continue; // because earlier loop may have undone the tag
6084  tile_ptr->tagged = false;
6085  for (TiledJet * jetI = tile_ptr->head; jetI != NULL; jetI = jetI->next) {
6086  if (jetI->NN == jetA || (jetI->NN == jetB && jetB != NULL)) {
6087  _set_NN(jetI, jets_for_minheap);
6088  }
6089  }
6090  }
6091  while (jets_for_minheap.size() > 0) {
6092  TiledJet * jetI = jets_for_minheap.back();
6093  jets_for_minheap.pop_back();
6094  minheap.update(jetI-head, _bj_diJ(jetI));
6095  jetI->label_minheap_update_done();
6096  Tile & tile_I = _tiles[jetI->tile_index];
6097  if (tile_I.max_NN_dist < jetI->NN_dist) tile_I.max_NN_dist = jetI->NN_dist;
6098  }
6099  n--;
6100  }
6101  delete[] briefjets;
6102 }
6104 #include <iomanip>
6105 #include <limits>
6106 #include <cmath>
6107 using namespace std;
6108 FJCORE_BEGIN_NAMESPACE // defined in fastjet/internal/base.hh
6109 TilingExtent::TilingExtent(ClusterSequence & cs) {
6110  _determine_rapidity_extent(cs.jets());
6111 }
6112 TilingExtent::TilingExtent(const vector<PseudoJet> &particles) {
6113  _determine_rapidity_extent(particles);
6114 }
6115 void TilingExtent::_determine_rapidity_extent(const vector<PseudoJet> & particles) {
6116  int nrap = 20;
6117  int nbins = 2*nrap;
6118  vector<double> counts(nbins, 0);
6119  _minrap = numeric_limits<double>::max();
6120  _maxrap = -numeric_limits<double>::max();
6121  int ibin;
6122  for (unsigned i = 0; i < particles.size(); i++) {
6123  if (particles[i].E() == abs(particles[i].pz())) continue;
6124  double rap = particles[i].rap();
6125  if (rap < _minrap) _minrap = rap;
6126  if (rap > _maxrap) _maxrap = rap;
6127  ibin = int(rap+nrap);
6128  if (ibin < 0) ibin = 0;
6129  if (ibin >= nbins) ibin = nbins - 1;
6130  counts[ibin]++;
6131  }
6132  double max_in_bin = 0;
6133  for (ibin = 0; ibin < nbins; ibin++) {
6134  if (max_in_bin < counts[ibin]) max_in_bin = counts[ibin];
6135  }
6136  const double allowed_max_fraction = 0.25;
6137  const double min_multiplicity = 4;
6138  double allowed_max_cumul = floor(max(max_in_bin * allowed_max_fraction, min_multiplicity));
6139  if (allowed_max_cumul > max_in_bin) allowed_max_cumul = max_in_bin;
6140  double cumul_lo = 0;
6141  _cumul2 = 0;
6142  for (ibin = 0; ibin < nbins; ibin++) {
6143  cumul_lo += counts[ibin];
6144  if (cumul_lo >= allowed_max_cumul) {
6145  double y = ibin-nrap;
6146  if (y > _minrap) _minrap = y;
6147  break;
6148  }
6149  }
6150  assert(ibin != nbins); // internal consistency check that you found a bin
6151  _cumul2 += cumul_lo*cumul_lo;
6152  int ibin_lo = ibin;
6153  double cumul_hi = 0;
6154  for (ibin = nbins-1; ibin >= 0; ibin--) {
6155  cumul_hi += counts[ibin];
6156  if (cumul_hi >= allowed_max_cumul) {
6157  double y = ibin-nrap+1; // +1 here is the rapidity bin width
6158  if (y < _maxrap) _maxrap = y;
6159  break;
6160  }
6161  }
6162  assert(ibin >= 0); // internal consistency check that you found a bin
6163  int ibin_hi = ibin;
6164  assert(ibin_hi >= ibin_lo);
6165  if (ibin_hi == ibin_lo) {
6166  _cumul2 = pow(double(cumul_lo + cumul_hi - counts[ibin_hi]), 2);
6167  } else {
6168  _cumul2 += cumul_hi*cumul_hi;
6169  for (ibin = ibin_lo+1; ibin < ibin_hi; ibin++) {
6170  _cumul2 += counts[ibin]*counts[ibin];
6171  }
6172  }
6173 }