Line data Source code
1 : // fjcore -- extracted from FastJet v3.0.5 (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 : // The results are expected to be identical to those obtained by linking to
13 : // the full FastJet distribution.
14 : //
15 : // NOTE THAT, IN ORDER TO MAKE IT POSSIBLE FOR FJCORE AND THE FULL FASTJET
16 : // TO COEXIST, THE FORMER USES THE "fjcore" NAMESPACE INSTEAD OF "fastjet".
17 : //
18 : // In particular, fjcore provides:
19 : //
20 : // - access to all native pp and ee algorithms, kt, anti-kt, C/A.
21 : // For C/A, the NlnN method is available, while anti-kt and kt
22 : // are limited to the N^2 one (still the fastest for N < 20k particles)
23 : // - access to selectors, for implementing cuts and selections
24 : // - access to all functionalities related to pseudojets (e.g. a jet's
25 : // structure or user-defined information)
26 : //
27 : // Instead, it does NOT provide:
28 : //
29 : // - jet areas functionality
30 : // - background estimation
31 : // - access to other algorithms via plugins
32 : // - interface to CGAL
33 : // - fastjet tools, e.g. filters, taggers
34 : //
35 : // If these functionalities are needed, the full FastJet installation must be
36 : // used. The code will be fully compatible, with the sole replacement of the
37 : // header files and of the fjcore namespace with the fastjet one.
38 : //
39 : // fjcore.hh and fjcore.cc are not meant to be human-readable.
40 : // For documentation, see the full FastJet manual and doxygen at http://fastjet.fr
41 : //
42 : // Like FastJet, fjcore is released under the terms of the GNU General Public
43 : // License version 2 (GPLv2). If you use this code as part of work towards a
44 : // scientific publication, whether directly or contained within another program
45 : // (e.g. Delphes, SpartyJet, Rivet, LHC collaboration software frameworks,
46 : // etc.), you should include a citation to
47 : //
48 : // EPJC72(2012)1896 [arXiv:1111.6097] (FastJet User Manual)
49 : // and, optionally, Phys.Lett.B641 (2006) 57 [arXiv:hep-ph/0512210]
50 : //
51 : // Copyright (c) 2005-2013, Matteo Cacciari, Gavin P. Salam and Gregory Soyez
52 : //
53 : //----------------------------------------------------------------------
54 : // This file is part of FastJet.
55 : //
56 : // FastJet is free software; you can redistribute it and/or modify
57 : // it under the terms of the GNU General Public License as published by
58 : // the Free Software Foundation; either version 2 of the License, or
59 : // (at your option) any later version.
60 : //
61 : // The algorithms that underlie FastJet have required considerable
62 : // development and are described in hep-ph/0512210. If you use
63 : // FastJet as part of work towards a scientific publication, please
64 : // include a citation to the FastJet paper.
65 : //
66 : // FastJet is distributed in the hope that it will be useful,
67 : // but WITHOUT ANY WARRANTY; without even the implied warranty of
68 : // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
69 : // GNU General Public License for more details.
70 : //
71 : // You should have received a copy of the GNU General Public License
72 : // along with FastJet. If not, see <http://www.gnu.org/licenses/>.
73 : //----------------------------------------------------------------------
74 : //
75 : //#include "fjcore.hh"
76 : // For inclusion in Pythia8 line above is replaced by line below.
77 : #include "Pythia8/FJcore.h"
78 : #ifndef __FJCORE_VERSION_HH__
79 : #define __FJCORE_VERSION_HH__
80 : #include<string>
81 : FJCORE_BEGIN_NAMESPACE
82 : const char* fastjet_version = FJCORE_PACKAGE_VERSION;
83 : FJCORE_END_NAMESPACE
84 : #endif // __FJCORE_VERSION_HH__
85 : #ifndef __FJCORE_CLUSTERQUENCE_N2_ICC__
86 : #define __FJCORE_CLUSTERQUENCE_N2_ICC__
87 : FJCORE_BEGIN_NAMESPACE // defined in fastjet/internal/base.hh
88 : template<class BJ> void ClusterSequence::_simple_N2_cluster() {
89 0 : int n = _jets.size();
90 0 : BJ * briefjets = new BJ[n];
91 0 : BJ * jetA = briefjets, * jetB;
92 0 : for (int i = 0; i< n; i++) {
93 0 : _bj_set_jetinfo(jetA, i);
94 0 : jetA++; // move on to next entry of briefjets
95 : }
96 0 : BJ * tail = jetA; // a semaphore for the end of briefjets
97 : BJ * head = briefjets; // a nicer way of naming start
98 0 : for (jetA = head + 1; jetA != tail; jetA++) {
99 0 : _bj_set_NN_crosscheck(jetA, head, jetA);
100 : }
101 0 : double * diJ = new double[n];
102 0 : jetA = head;
103 0 : for (int i = 0; i < n; i++) {
104 0 : diJ[i] = _bj_diJ(jetA);
105 0 : jetA++; // have jetA follow i
106 : }
107 0 : int history_location = n-1;
108 0 : while (tail != head) {
109 0 : double diJ_min = diJ[0];
110 : int diJ_min_jet = 0;
111 0 : for (int i = 1; i < n; i++) {
112 0 : if (diJ[i] < diJ_min) {diJ_min_jet = i; diJ_min = diJ[i];}
113 : }
114 0 : history_location++;
115 0 : jetA = & briefjets[diJ_min_jet];
116 0 : jetB = static_cast<BJ *>(jetA->NN);
117 0 : diJ_min *= _invR2;
118 0 : if (jetB != NULL) {
119 0 : if (jetA < jetB) {std::swap(jetA,jetB);}
120 0 : int nn; // new jet index
121 0 : _do_ij_recombination_step(jetA->_jets_index, jetB->_jets_index, diJ_min, nn);
122 0 : _bj_set_jetinfo(jetB, nn);
123 0 : } else {
124 0 : _do_iB_recombination_step(jetA->_jets_index, diJ_min);
125 : }
126 0 : tail--; n--;
127 0 : *jetA = *tail;
128 0 : diJ[jetA - head] = diJ[tail-head];
129 0 : for (BJ * jetI = head; jetI != tail; jetI++) {
130 0 : if (jetI->NN == jetA || jetI->NN == jetB) {
131 0 : _bj_set_NN_nocross(jetI, head, tail);
132 0 : diJ[jetI-head] = _bj_diJ(jetI); // update diJ
133 0 : }
134 0 : if (jetB != NULL) {
135 0 : double dist = _bj_dist(jetI,jetB);
136 0 : if (dist < jetI->NN_dist) {
137 0 : if (jetI != jetB) {
138 0 : jetI->NN_dist = dist;
139 0 : jetI->NN = jetB;
140 0 : diJ[jetI-head] = _bj_diJ(jetI); // update diJ...
141 0 : }
142 : }
143 0 : if (dist < jetB->NN_dist) {
144 0 : if (jetI != jetB) {
145 0 : jetB->NN_dist = dist;
146 0 : jetB->NN = jetI;}
147 : }
148 0 : }
149 0 : if (jetI->NN == tail) {jetI->NN = jetA;}
150 : }
151 0 : if (jetB != NULL) {diJ[jetB-head] = _bj_diJ(jetB);}
152 0 : }
153 0 : delete[] diJ;
154 0 : delete[] briefjets;
155 0 : }
156 : FJCORE_END_NAMESPACE
157 : #endif // __FJCORE_CLUSTERQUENCE_N2_ICC__
158 : #ifndef __FJCORE_DYNAMICNEARESTNEIGHBOURS_HH__
159 : #define __FJCORE_DYNAMICNEARESTNEIGHBOURS_HH__
160 : #include<vector>
161 : #include<string>
162 : #include<iostream>
163 : #include<sstream>
164 : #include<cassert>
165 : FJCORE_BEGIN_NAMESPACE // defined in fastjet/internal/base.hh
166 : class EtaPhi {
167 : public:
168 : double first, second;
169 0 : EtaPhi() {}
170 0 : EtaPhi(double a, double b) {first = a; second = b;}
171 : void sanitize() {
172 0 : if (second < 0) second += twopi;
173 0 : if (second >= twopi) second -= twopi;
174 0 : }
175 : };
176 : class DnnError {
177 : public:
178 : DnnError() {;};
179 : DnnError(const std::string & message_in) {
180 : _message = message_in; std::cerr << message_in << std::endl;};
181 : std::string message() const {return _message;};
182 : private:
183 : std::string _message;
184 : };
185 : class DynamicNearestNeighbours {
186 : public:
187 : virtual int NearestNeighbourIndex(const int & ii) const = 0;
188 : virtual double NearestNeighbourDistance(const int & ii) const = 0;
189 : virtual bool Valid(const int & index) const = 0;
190 : virtual void RemoveAndAddPoints(const std::vector<int> & indices_to_remove,
191 : const std::vector<EtaPhi> & points_to_add,
192 : std::vector<int> & indices_added,
193 : std::vector<int> & indices_of_updated_neighbours) = 0;
194 : inline void RemovePoint (const int & index,
195 : std::vector<int> & indices_of_updated_neighbours) {
196 : std::vector<int> indices_added;
197 : std::vector<EtaPhi> points_to_add;
198 : std::vector<int> indices_to_remove(1);
199 : indices_to_remove[0] = index;
200 : RemoveAndAddPoints(indices_to_remove, points_to_add, indices_added,
201 : indices_of_updated_neighbours
202 : );};
203 : inline void RemoveCombinedAddCombination(
204 : const int & index1, const int & index2,
205 : const EtaPhi & newpoint,
206 : int & index3,
207 : std::vector<int> & indices_of_updated_neighbours) {
208 : std::vector<int> indices_added(1);
209 : std::vector<EtaPhi> points_to_add(1);
210 : std::vector<int> indices_to_remove(2);
211 : indices_to_remove[0] = index1;
212 : indices_to_remove[1] = index2;
213 : points_to_add[0] = newpoint;
214 : RemoveAndAddPoints(indices_to_remove, points_to_add, indices_added,
215 : indices_of_updated_neighbours
216 : );
217 : index3 = indices_added[0];
218 : };
219 : virtual ~DynamicNearestNeighbours () {}
220 : };
221 : FJCORE_END_NAMESPACE
222 : #endif // __FJCORE_DYNAMICNEARESTNEIGHBOURS_HH__
223 : #ifndef __FJCORE_SEARCHTREE_HH__
224 : #define __FJCORE_SEARCHTREE_HH__
225 : #include<vector>
226 : #include<cassert>
227 : #include<cstddef>
228 : FJCORE_BEGIN_NAMESPACE // defined in fastjet/internal/base.hh
229 0 : template<class T> class SearchTree {
230 : public:
231 : class Node;
232 : class circulator;
233 : class const_circulator;
234 : SearchTree(const std::vector<T> & init);
235 0 : SearchTree(const std::vector<T> & init, unsigned int max_size);
236 : void remove(unsigned node_index);
237 : void remove(typename SearchTree::Node * node);
238 : void remove(typename SearchTree::circulator & circ);
239 : circulator insert(const T & value);
240 : const Node & operator[](int i) const {return _nodes[i];};
241 0 : unsigned int size() const {return _nodes.size() - _available_nodes.size();}
242 : void verify_structure();
243 : void verify_structure_linear() const;
244 : void verify_structure_recursive(const Node * , const Node * , const Node * ) const;
245 : void print_elements();
246 : #ifdef __FJCORE_SEARCHTREE_TRACK_DEPTH
247 : inline unsigned int max_depth() const {return _max_depth;};
248 : #else
249 : inline unsigned int max_depth() const {return 0;};
250 : #endif
251 : int loc(const Node * node) const ;
252 : Node * _find_predecessor(const Node *);
253 : Node * _find_successor(const Node *);
254 : const Node & operator[](unsigned int i) const {return _nodes[i];};
255 : const_circulator somewhere() const;
256 : circulator somewhere();
257 : private:
258 : void _initialize(const std::vector<T> & init);
259 : std::vector<Node> _nodes;
260 : std::vector<Node *> _available_nodes;
261 : Node * _top_node;
262 : unsigned int _n_removes;
263 : void _do_initial_connections(unsigned int this_one, unsigned int scale,
264 : unsigned int left_edge, unsigned int right_edge,
265 : unsigned int depth);
266 : #ifdef __FJCORE_SEARCHTREE_TRACK_DEPTH
267 : unsigned int _max_depth;
268 : #endif
269 : };
270 : template<class T> class SearchTree<T>::Node{
271 : public:
272 0 : Node() {}; /// default constructor
273 : bool treelinks_null() const {
274 0 : return ((parent==0) && (left==0) && (right==0));};
275 : inline void nullify_treelinks() {
276 0 : parent = NULL;
277 0 : left = NULL;
278 0 : right = NULL;
279 0 : };
280 : void reset_parents_link_to_me(Node * XX);
281 : T value;
282 : Node * left;
283 : Node * right;
284 : Node * parent;
285 : Node * successor;
286 : Node * predecessor;
287 : };
288 : template<class T> void SearchTree<T>::Node::reset_parents_link_to_me(typename SearchTree<T>::Node * XX) {
289 0 : if (parent == NULL) {return;}
290 0 : if (parent->right == this) {parent->right = XX;}
291 0 : else {parent->left = XX;}
292 0 : }
293 : template<class T> class SearchTree<T>::circulator{
294 : public:
295 : // Next line commented out, by author agreement, to avoid compiler warning.
296 : //template<class U> friend class SearchTree<U>::const_circulator;
297 : friend class SearchTree<T>;
298 0 : circulator() : _node(NULL) {}
299 0 : circulator(Node * node) : _node(node) {}
300 : const T * operator->() const {return &(_node->value);}
301 0 : T * operator->() {return &(_node->value);}
302 : const T & operator*() const {return _node->value;}
303 : T & operator*() {return _node->value;}
304 : circulator & operator++() {
305 0 : _node = _node->successor;
306 0 : return *this;}
307 : circulator operator++(int) {
308 0 : circulator tmp = *this;
309 0 : _node = _node->successor;
310 0 : return tmp;}
311 : circulator & operator--() {
312 : _node = _node->predecessor;
313 : return *this;}
314 : circulator operator--(int) {
315 0 : circulator tmp = *this;
316 0 : _node = _node->predecessor;
317 0 : return tmp;}
318 : circulator next() const {
319 0 : return circulator(_node->successor);}
320 : circulator previous() const {
321 : return circulator(_node->predecessor);}
322 0 : bool operator!=(const circulator & other) const {return other._node != _node;}
323 : bool operator==(const circulator & other) const {return other._node == _node;}
324 : private:
325 : Node * _node;
326 : };
327 : template<class T> class SearchTree<T>::const_circulator{
328 : public:
329 : const_circulator() : _node(NULL) {}
330 : const_circulator(const Node * node) : _node(node) {}
331 : const_circulator(const circulator & circ) :_node(circ._node) {}
332 : const T * operator->() {return &(_node->value);}
333 : const T & operator*() const {return _node->value;}
334 : const_circulator & operator++() {
335 : _node = _node->successor;
336 : return *this;}
337 : const_circulator operator++(int) {
338 : const_circulator tmp = *this;
339 : _node = _node->successor;
340 : return tmp;}
341 : const_circulator & operator--() {
342 : _node = _node->predecessor;
343 : return *this;}
344 : const_circulator operator--(int) {
345 : const_circulator tmp = *this;
346 : _node = _node->predecessor;
347 : return tmp;}
348 : const_circulator next() const {
349 : return const_circulator(_node->successor);}
350 : const_circulator previous() const {
351 : return const_circulator(_node->predecessor);}
352 : bool operator!=(const const_circulator & other) const {return other._node != _node;}
353 : bool operator==(const const_circulator & other) const {return other._node == _node;}
354 : private:
355 : const Node * _node;
356 : };
357 : template<class T> SearchTree<T>::SearchTree(const std::vector<T> & init,
358 : unsigned int max_size) :
359 0 : _nodes(max_size) {
360 0 : _available_nodes.reserve(max_size);
361 0 : _available_nodes.resize(max_size - init.size());
362 0 : for (unsigned int i = init.size(); i < max_size; i++) {
363 0 : _available_nodes[i-init.size()] = &(_nodes[i]);
364 : }
365 0 : _initialize(init);
366 0 : }
367 : template<class T> SearchTree<T>::SearchTree(const std::vector<T> & init) :
368 : _nodes(init.size()), _available_nodes(0) {
369 : _available_nodes.reserve(init.size());
370 : _initialize(init);
371 : }
372 : template<class T> void SearchTree<T>::_initialize(const std::vector<T> & init) {
373 0 : _n_removes = 0;
374 0 : unsigned n = init.size();
375 0 : assert(n>=1);
376 : #ifdef __FJCORE_SEARCHTREE_TRACK_DEPTH
377 : _max_depth = 0;
378 : #endif
379 0 : for (unsigned int i = 1; i<n; i++) {
380 0 : assert(!(init[i] < init[i-1]));
381 : }
382 0 : for(unsigned int i = 0; i < n; i++) {
383 0 : _nodes[i].value = init[i];
384 0 : _nodes[i].predecessor = (& (_nodes[i])) - 1;
385 0 : _nodes[i].successor = (& (_nodes[i])) + 1;
386 0 : _nodes[i].nullify_treelinks();
387 : }
388 0 : _nodes[0].predecessor = (& (_nodes[n-1]));
389 0 : _nodes[n-1].successor = (& (_nodes[0]));
390 0 : unsigned int scale = (n+1)/2;
391 0 : unsigned int top = std::min(n-1,scale);
392 0 : _nodes[top].parent = NULL;
393 0 : _top_node = &(_nodes[top]);
394 0 : _do_initial_connections(top, scale, 0, n, 0);
395 0 : }
396 : template<class T> inline int SearchTree<T>::loc(const Node * node) const {return node == NULL?
397 : -999 : node - &(_nodes[0]);}
398 : template<class T> void SearchTree<T>::_do_initial_connections(
399 : unsigned int this_one,
400 : unsigned int scale,
401 : unsigned int left_edge,
402 : unsigned int right_edge,
403 : unsigned int depth
404 : ) {
405 : #ifdef __FJCORE_SEARCHTREE_TRACK_DEPTH
406 : _max_depth = max(depth, _max_depth);
407 : #endif
408 0 : unsigned int ref_new_scale = (scale+1)/2;
409 : unsigned new_scale = ref_new_scale;
410 : bool did_child = false;
411 0 : while(true) {
412 0 : int left = this_one - new_scale; // be careful here to use signed int...
413 0 : if (left >= static_cast<int>(left_edge)
414 0 : && _nodes[left].treelinks_null() ) {
415 0 : _nodes[left].parent = &(_nodes[this_one]);
416 0 : _nodes[this_one].left = &(_nodes[left]);
417 0 : _do_initial_connections(left, new_scale, left_edge, this_one, depth+1);
418 : did_child = true;
419 0 : break;
420 : }
421 : unsigned int old_new_scale = new_scale;
422 0 : new_scale = (old_new_scale + 1)/2;
423 0 : if (new_scale == old_new_scale) break;
424 0 : }
425 0 : if (!did_child) {_nodes[this_one].left = NULL;}
426 : new_scale = ref_new_scale;
427 : did_child = false;
428 0 : while(true) {
429 0 : unsigned int right = this_one + new_scale;
430 0 : if (right < right_edge && _nodes[right].treelinks_null()) {
431 0 : _nodes[right].parent = &(_nodes[this_one]);
432 0 : _nodes[this_one].right = &(_nodes[right]);
433 0 : _do_initial_connections(right, new_scale, this_one+1,right_edge,depth+1);
434 : did_child = true;
435 0 : break;
436 : }
437 : unsigned int old_new_scale = new_scale;
438 0 : new_scale = (old_new_scale + 1)/2;
439 0 : if (new_scale == old_new_scale) break;
440 0 : }
441 0 : if (!did_child) {_nodes[this_one].right = NULL;}
442 0 : }
443 : template<class T> void SearchTree<T>::remove(unsigned int node_index) {
444 : remove(&(_nodes[node_index]));
445 : }
446 : template<class T> void SearchTree<T>::remove(circulator & circ) {
447 0 : remove(circ._node);
448 0 : }
449 : template<class T> void SearchTree<T>::remove(typename SearchTree<T>::Node * node) {
450 0 : assert(size() > 1); // switch this to throw...?
451 0 : assert(!node->treelinks_null());
452 0 : node->predecessor->successor = node->successor;
453 0 : node->successor->predecessor = node->predecessor;
454 0 : if (node->left == NULL && node->right == NULL) {
455 0 : node->reset_parents_link_to_me(NULL);
456 0 : } else if (node->left != NULL && node->right == NULL){
457 0 : node->reset_parents_link_to_me(node->left);
458 0 : node->left->parent = node->parent;
459 0 : if (_top_node == node) {_top_node = node->left;}
460 0 : } else if (node->left == NULL && node->right != NULL){
461 0 : node->reset_parents_link_to_me(node->right);
462 0 : node->right->parent = node->parent;
463 0 : if (_top_node == node) {_top_node = node->right;}
464 : } else {
465 : Node * replacement;
466 0 : bool use_predecessor = (_n_removes % 2 == 1);
467 0 : if (use_predecessor) {
468 0 : replacement = node->predecessor;
469 0 : assert(replacement->right == NULL); // guaranteed if it's our predecessor
470 0 : if (replacement != node->left) {
471 0 : if (replacement->left != NULL) {
472 0 : replacement->left->parent = replacement->parent;}
473 0 : replacement->reset_parents_link_to_me(replacement->left);
474 0 : replacement->left = node->left;
475 0 : }
476 0 : replacement->parent = node->parent;
477 0 : replacement->right = node->right;
478 0 : } else {
479 0 : replacement = node->successor;
480 0 : assert(replacement->left == NULL); // guaranteed if it's our successor
481 0 : if (replacement != node->right) {
482 0 : if (replacement->right != NULL) {
483 0 : replacement->right->parent = replacement->parent;}
484 0 : replacement->reset_parents_link_to_me(replacement->right);
485 0 : replacement->right = node->right;
486 0 : }
487 0 : replacement->parent = node->parent;
488 0 : replacement->left = node->left;
489 : }
490 0 : node->reset_parents_link_to_me(replacement);
491 0 : if (node->left != replacement) {node->left->parent = replacement;}
492 0 : if (node->right != replacement) {node->right->parent = replacement;}
493 0 : if (_top_node == node) {_top_node = replacement;}
494 : }
495 0 : node->nullify_treelinks();
496 0 : node->predecessor = NULL;
497 0 : node->successor = NULL;
498 0 : _n_removes++;
499 0 : _available_nodes.push_back(node);
500 0 : }
501 : template<class T> typename SearchTree<T>::circulator SearchTree<T>::insert(const T & value) {
502 0 : assert(_available_nodes.size() > 0);
503 0 : Node * node = _available_nodes.back();
504 0 : _available_nodes.pop_back();
505 0 : node->value = value;
506 0 : Node * location = _top_node;
507 : Node * old_location = NULL;
508 : bool on_left = true; // (init not needed -- but soothes g++4)
509 : #ifdef __FJCORE_SEARCHTREE_TRACK_DEPTH
510 : unsigned int depth = 0;
511 : #endif
512 0 : while(location != NULL) {
513 : #ifdef __FJCORE_SEARCHTREE_TRACK_DEPTH
514 : depth++;
515 : #endif
516 : old_location = location;
517 0 : on_left = value < location->value;
518 0 : if (on_left) {location = location->left;}
519 0 : else {location = location->right;}
520 : }
521 : #ifdef __FJCORE_SEARCHTREE_TRACK_DEPTH
522 : _max_depth = max(depth, _max_depth);
523 : #endif
524 0 : node->parent = old_location;
525 0 : if (on_left) {node->parent->left = node;}
526 0 : else {node->parent->right = node;}
527 0 : node->left = NULL;
528 0 : node->right = NULL;
529 0 : node->predecessor = _find_predecessor(node);
530 0 : if (node->predecessor != NULL) {
531 0 : node->successor = node->predecessor->successor;
532 0 : node->predecessor->successor = node;
533 0 : node->successor->predecessor = node;
534 0 : } else {
535 0 : node->successor = _find_successor(node);
536 0 : assert(node->successor != NULL); // can only happen if we're sole element
537 0 : node->predecessor = node->successor->predecessor;
538 0 : node->successor->predecessor = node;
539 0 : node->predecessor->successor = node;
540 : }
541 0 : return circulator(node);
542 0 : }
543 : template<class T> void SearchTree<T>::verify_structure() {
544 : verify_structure_linear();
545 : const Node * left_limit = _top_node;
546 : while (left_limit->left != NULL) {left_limit = left_limit->left;}
547 : const Node * right_limit = _top_node;
548 : while (right_limit->right != NULL) {right_limit = right_limit->right;}
549 : verify_structure_recursive(_top_node, left_limit, right_limit);
550 : }
551 : template<class T> void SearchTree<T>::verify_structure_recursive(
552 : const typename SearchTree<T>::Node * element,
553 : const typename SearchTree<T>::Node * left_limit,
554 : const typename SearchTree<T>::Node * right_limit) const {
555 : assert(!(element->value < left_limit->value));
556 : assert(!(right_limit->value < element->value));
557 : const Node * left = element->left;
558 : if (left != NULL) {
559 : assert(!(element->value < left->value));
560 : if (left != left_limit) {
561 : verify_structure_recursive(left, left_limit, element);}
562 : }
563 : const Node * right = element->right;
564 : if (right != NULL) {
565 : assert(!(right->value < element->value));
566 : if (right != right_limit) {
567 : verify_structure_recursive(right, element, right_limit);}
568 : }
569 : }
570 : template<class T> void SearchTree<T>::verify_structure_linear() const {
571 : unsigned n_top = 0;
572 : unsigned n_null = 0;
573 : for(unsigned i = 0; i < _nodes.size(); i++) {
574 : const typename SearchTree<T>::Node * node = &(_nodes[i]);
575 : if (node->treelinks_null()) {n_null++; continue;}
576 : if (node->parent == NULL) {
577 : n_top++;
578 : } else {
579 : assert((node->parent->left == node) ^ (node->parent->right == node));
580 : }
581 : if (node->left != NULL) {
582 : assert(!(node->value < node->left->value ));}
583 : if (node->right != NULL) {
584 : assert(!(node->right->value < node->value ));}
585 : }
586 : assert(n_top == 1 || (n_top == 0 && size() <= 1) );
587 : assert(n_null == _available_nodes.size() ||
588 : (n_null == _available_nodes.size() + 1 && size() == 1));
589 : }
590 : template<class T> typename SearchTree<T>::Node * SearchTree<T>::_find_predecessor(const typename SearchTree<T>::Node * node) {
591 : typename SearchTree<T>::Node * newnode;
592 0 : if (node->left != NULL) {
593 : newnode = node->left;
594 0 : while(newnode->right != NULL) {newnode = newnode->right;}
595 0 : return newnode;
596 : } else {
597 : const typename SearchTree<T>::Node * lastnode = node;
598 0 : newnode = node->parent;
599 0 : while(newnode != NULL) {
600 0 : if (newnode->right == lastnode) {return newnode;}
601 : lastnode = newnode;
602 0 : newnode = newnode->parent;
603 : }
604 0 : return newnode;
605 : }
606 0 : }
607 : template<class T> typename SearchTree<T>::Node * SearchTree<T>::_find_successor(const typename SearchTree<T>::Node * node) {
608 : typename SearchTree<T>::Node * newnode;
609 0 : if (node->right != NULL) {
610 : newnode = node->right;
611 0 : while(newnode->left != NULL) {newnode = newnode->left;}
612 0 : return newnode;
613 : } else {
614 : const typename SearchTree<T>::Node * lastnode = node;
615 0 : newnode = node->parent;
616 0 : while(newnode != NULL) {
617 0 : if (newnode->left == lastnode) {return newnode;}
618 : lastnode = newnode;
619 0 : newnode = newnode->parent;
620 : }
621 0 : return newnode;
622 : }
623 0 : }
624 : template<class T> void SearchTree<T>::print_elements() {
625 : typename SearchTree<T>::Node * base_node = &(_nodes[0]);
626 : typename SearchTree<T>::Node * node = base_node;
627 : int n = _nodes.size();
628 : for(; node - base_node < n ; node++) {
629 : 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);
630 : }
631 : }
632 : template<class T> typename SearchTree<T>::circulator SearchTree<T>::somewhere() {
633 0 : return circulator(_top_node);
634 : }
635 : template<class T> typename SearchTree<T>::const_circulator SearchTree<T>::somewhere() const {
636 : return const_circulator(_top_node);
637 : }
638 : FJCORE_END_NAMESPACE
639 : #endif // __FJCORE_SEARCHTREE_HH__
640 : #ifndef __FJCORE_MINHEAP__HH__
641 : #define __FJCORE_MINHEAP__HH__
642 : #include<vector>
643 : #include<cassert>
644 : #include<memory>
645 : #include<limits>
646 : FJCORE_BEGIN_NAMESPACE // defined in fastjet/internal/base.hh
647 0 : class MinHeap {
648 : public:
649 : MinHeap (const std::vector<double> & values, unsigned int max_size) :
650 0 : _heap(max_size) {_initialise(values);};
651 : MinHeap (const std::vector<double> & values) :
652 0 : _heap(values.size()) {_initialise(values);};
653 : inline unsigned int minloc() const {
654 0 : return (_heap[0].minloc) - &(_heap[0]);};
655 0 : inline double minval() const {return _heap[0].minloc->value;};
656 : inline double operator[](int i) const {return _heap[i].value;};
657 : void remove(unsigned int loc) {
658 0 : update(loc,std::numeric_limits<double>::max());};
659 : void update(unsigned int, double);
660 : private:
661 : struct ValueLoc{
662 : double value;
663 : ValueLoc * minloc;
664 : };
665 : std::vector<ValueLoc> _heap;
666 : void _initialise(const std::vector<double> & values);
667 : };
668 : FJCORE_END_NAMESPACE
669 : #endif // __FJCORE_MINHEAP__HH__
670 : #ifndef __FJCORE_CLOSESTPAIR2DBASE__HH__
671 : #define __FJCORE_CLOSESTPAIR2DBASE__HH__
672 : #include<vector>
673 : FJCORE_BEGIN_NAMESPACE // defined in fastjet/internal/base.hh
674 : class Coord2D {
675 : public:
676 : double x, y;
677 0 : Coord2D() {};
678 0 : Coord2D(double a, double b): x(a), y(b) {};
679 : Coord2D operator-(const Coord2D & other) const {
680 0 : return Coord2D(x - other.x, y - other.y);};
681 : Coord2D operator+(const Coord2D & other) const {
682 : return Coord2D(x + other.x, y + other.y);};
683 : Coord2D operator*(double factor) const {return Coord2D(factor*x,factor*y);};
684 : friend Coord2D operator*(double factor, const Coord2D & coord) {
685 : return Coord2D(factor*coord.x,factor*coord.y);
686 : }
687 : Coord2D operator/(double divisor) const {
688 0 : return Coord2D(x / divisor, y / divisor);};
689 : friend double distance2(const Coord2D & a, const Coord2D & b) {
690 : double dx = a.x - b.x, dy = a.y-b.y;
691 : return dx*dx+dy*dy;
692 : };
693 : double distance2(const Coord2D & b) const {
694 0 : double dx = x - b.x, dy = y-b.y;
695 0 : return dx*dx+dy*dy;
696 : };
697 : };
698 0 : class ClosestPair2DBase {
699 : public:
700 : virtual void closest_pair(unsigned int & ID1, unsigned int & ID2,
701 : double & distance2) const = 0;
702 : virtual void remove(unsigned int ID) = 0;
703 : virtual unsigned int insert(const Coord2D & position) = 0;
704 : virtual unsigned int replace(unsigned int ID1, unsigned int ID2,
705 : const Coord2D & position) {
706 0 : remove(ID1);
707 0 : remove(ID2);
708 0 : unsigned new_ID = insert(position);
709 0 : return(new_ID);
710 : };
711 : virtual void replace_many(const std::vector<unsigned int> & IDs_to_remove,
712 : const std::vector<Coord2D> & new_positions,
713 : std::vector<unsigned int> & new_IDs) {
714 0 : for(unsigned i = 0; i < IDs_to_remove.size(); i++) {
715 0 : remove(IDs_to_remove[i]);}
716 0 : new_IDs.resize(0);
717 0 : for(unsigned i = 0; i < new_positions.size(); i++) {
718 0 : new_IDs.push_back(insert(new_positions[i]));}
719 0 : }
720 : virtual unsigned int size() = 0;
721 0 : virtual ~ClosestPair2DBase() {};
722 : };
723 : FJCORE_END_NAMESPACE
724 : #endif // __FJCORE_CLOSESTPAIR2DBASE__HH__
725 : #ifndef __FJCORE_CLOSESTPAIR2D__HH__
726 : #define __FJCORE_CLOSESTPAIR2D__HH__
727 : #include<vector>
728 : #include<stack>
729 : #include<iostream>
730 : FJCORE_BEGIN_NAMESPACE // defined in fastjet/internal/base.hh
731 0 : class ClosestPair2D : public ClosestPair2DBase {
732 : public:
733 0 : ClosestPair2D(const std::vector<Coord2D> & positions,
734 0 : const Coord2D & left_corner, const Coord2D & right_corner) {
735 0 : _initialize(positions, left_corner, right_corner, positions.size());
736 0 : };
737 : ClosestPair2D(const std::vector<Coord2D> & positions,
738 : const Coord2D & left_corner, const Coord2D & right_corner,
739 : const unsigned int max_size) {
740 : _initialize(positions, left_corner, right_corner, max_size);
741 : };
742 : void closest_pair(unsigned int & ID1, unsigned int & ID2,
743 : double & distance2) const;
744 : void remove(unsigned int ID);
745 : unsigned int insert(const Coord2D &);
746 : virtual unsigned int replace(unsigned int ID1, unsigned int ID2,
747 : const Coord2D & position);
748 : virtual void replace_many(const std::vector<unsigned int> & IDs_to_remove,
749 : const std::vector<Coord2D> & new_positions,
750 : std::vector<unsigned int> & new_IDs);
751 : inline void print_tree_depths(std::ostream & outdev) const {
752 : outdev << _trees[0]->max_depth() << " "
753 : << _trees[1]->max_depth() << " "
754 : << _trees[2]->max_depth() << "\n";
755 : };
756 : unsigned int size();
757 : private:
758 : void _initialize(const std::vector<Coord2D> & positions,
759 : const Coord2D & left_corner, const Coord2D & right_corner,
760 : const unsigned int max_size);
761 : static const unsigned int _nshift = 3;
762 : class Point; // will be defined below
763 0 : template<class T> class triplet {
764 : public:
765 : inline const T & operator[](unsigned int i) const {return _contents[i];};
766 0 : inline T & operator[](unsigned int i) {return _contents[i];};
767 : private:
768 : T _contents[_nshift];
769 : };
770 : class Shuffle {
771 : public:
772 : unsigned int x, y;
773 : Point * point;
774 : bool operator<(const Shuffle &) const;
775 0 : void operator+=(unsigned int shift) {x += shift; y+= shift;};
776 : };
777 : typedef SearchTree<Shuffle> Tree;
778 : typedef Tree::circulator circulator;
779 : typedef Tree::const_circulator const_circulator;
780 : triplet<std::auto_ptr<Tree> > _trees;
781 : std::auto_ptr<MinHeap> _heap;
782 : std::vector<Point> _points;
783 : std::stack<Point *> _available_points;
784 : std::vector<Point *> _points_under_review;
785 : static const unsigned int _remove_heap_entry = 1;
786 : static const unsigned int _review_heap_entry = 2;
787 : static const unsigned int _review_neighbour = 4;
788 : void _add_label(Point * point, unsigned int review_flag);
789 : void _set_label(Point * point, unsigned int review_flag);
790 : void _deal_with_points_to_review();
791 : void _remove_from_search_tree(Point * point_to_remove);
792 : void _insert_into_search_tree(Point * new_point);
793 : void _point2shuffle(Point & , Shuffle & , unsigned int shift);
794 : Coord2D _left_corner;
795 : double _range;
796 : int _ID(const Point *) const;
797 : triplet<unsigned int> _shifts; // absolute shifts
798 : triplet<unsigned int> _rel_shifts; // shifts relative to previous shift
799 : unsigned int _cp_search_range;
800 : };
801 0 : class ClosestPair2D::Point {
802 : public:
803 : Coord2D coord;
804 : Point * neighbour;
805 : double neighbour_dist2;
806 : triplet<circulator> circ;
807 : unsigned int review_flag;
808 : double distance2(const Point & other) const {
809 0 : return coord.distance2(other.coord);
810 : };
811 : };
812 : inline bool floor_ln2_less(unsigned x, unsigned y) {
813 0 : if (x>y) return false;
814 0 : return (x < (x^y)); // beware of operator precedence...
815 0 : }
816 : inline int ClosestPair2D::_ID(const Point * point) const {
817 0 : return point - &(_points[0]);
818 : }
819 : inline unsigned int ClosestPair2D::size() {
820 0 : return _points.size() - _available_points.size();
821 : }
822 : FJCORE_END_NAMESPACE
823 : #endif // __FJCORE_CLOSESTPAIR2D__HH__
824 : #include<limits>
825 : #include<iostream>
826 : #include<iomanip>
827 : #include<algorithm>
828 : FJCORE_BEGIN_NAMESPACE // defined in fastjet/internal/base.hh
829 : // Next line commented out, by author agreement, to avoid compiler warning.
830 : //const unsigned int huge_unsigned = 4294967295U;
831 : const unsigned int twopow31 = 2147483648U;
832 : using namespace std;
833 : void ClosestPair2D::_point2shuffle(Point & point, Shuffle & shuffle,
834 : unsigned int shift) {
835 0 : Coord2D renorm_point = (point.coord - _left_corner)/_range;
836 0 : assert(renorm_point.x >=0);
837 0 : assert(renorm_point.x <=1);
838 0 : assert(renorm_point.y >=0);
839 0 : assert(renorm_point.y <=1);
840 0 : shuffle.x = static_cast<unsigned int>(twopow31 * renorm_point.x) + shift;
841 0 : shuffle.y = static_cast<unsigned int>(twopow31 * renorm_point.y) + shift;
842 0 : shuffle.point = &point;
843 0 : }
844 : bool ClosestPair2D::Shuffle::operator<(const Shuffle & q) const {
845 0 : if (floor_ln2_less(x ^ q.x, y ^ q.y)) {
846 0 : return (y < q.y);
847 : } else {
848 0 : return (x < q.x);
849 : }
850 0 : }
851 : void ClosestPair2D::_initialize(const std::vector<Coord2D> & positions,
852 : const Coord2D & left_corner,
853 : const Coord2D & right_corner,
854 : unsigned int max_size) {
855 0 : unsigned int n_positions = positions.size();
856 0 : assert(max_size >= n_positions);
857 0 : _points.resize(max_size);
858 0 : for (unsigned int i = n_positions; i < max_size; i++) {
859 0 : _available_points.push(&(_points[i]));
860 : }
861 0 : _left_corner = left_corner;
862 0 : _range = max((right_corner.x - left_corner.x),
863 0 : (right_corner.y - left_corner.y));
864 0 : vector<Shuffle> shuffles(n_positions);
865 0 : for (unsigned int i = 0; i < n_positions; i++) {
866 0 : _points[i].coord = positions[i];
867 0 : _points[i].neighbour_dist2 = numeric_limits<double>::max();
868 0 : _points[i].review_flag = 0;
869 0 : _point2shuffle(_points[i], shuffles[i], 0);
870 : }
871 0 : for (unsigned ishift = 0; ishift < _nshift; ishift++) {
872 0 : _shifts[ishift] = static_cast<unsigned int>(((twopow31*1.0)*ishift)/_nshift);
873 0 : if (ishift == 0) {_rel_shifts[ishift] = 0;}
874 0 : else {_rel_shifts[ishift] = _shifts[ishift] - _shifts[ishift-1];}
875 : }
876 0 : _cp_search_range = 30;
877 0 : _points_under_review.reserve(_nshift * _cp_search_range);
878 0 : for (unsigned int ishift = 0; ishift < _nshift; ishift++) {
879 0 : if (ishift > 0) {
880 0 : unsigned rel_shift = _rel_shifts[ishift];
881 0 : for (unsigned int i = 0; i < shuffles.size(); i++) {
882 0 : shuffles[i] += rel_shift; }
883 0 : }
884 0 : sort(shuffles.begin(), shuffles.end());
885 0 : _trees[ishift] = auto_ptr<Tree>(new Tree(shuffles, max_size));
886 0 : circulator circ = _trees[ishift]->somewhere(), start=circ;
887 0 : unsigned int CP_range = min(_cp_search_range, n_positions-1);
888 0 : do {
889 0 : Point * this_point = circ->point;
890 0 : this_point->circ[ishift] = circ;
891 0 : circulator other = circ;
892 0 : for (unsigned i=0; i < CP_range; i++) {
893 0 : ++other;
894 0 : double dist2 = this_point->distance2(*other->point);
895 0 : if (dist2 < this_point->neighbour_dist2) {
896 0 : this_point->neighbour_dist2 = dist2;
897 0 : this_point->neighbour = other->point;
898 0 : }
899 : }
900 0 : } while (++circ != start);
901 0 : }
902 0 : vector<double> mindists2(n_positions);
903 0 : for (unsigned int i = 0; i < n_positions; i++) {
904 0 : mindists2[i] = _points[i].neighbour_dist2;}
905 0 : _heap = auto_ptr<MinHeap>(new MinHeap(mindists2, max_size));
906 0 : }
907 : void ClosestPair2D::closest_pair(unsigned int & ID1, unsigned int & ID2,
908 : double & distance2) const {
909 0 : ID1 = _heap->minloc();
910 0 : ID2 = _ID(_points[ID1].neighbour);
911 0 : distance2 = _points[ID1].neighbour_dist2;
912 0 : if (ID1 > ID2) std::swap(ID1,ID2);
913 0 : }
914 : inline void ClosestPair2D::_add_label(Point * point, unsigned int review_flag) {
915 0 : if (point->review_flag == 0) _points_under_review.push_back(point);
916 0 : point->review_flag |= review_flag;
917 0 : }
918 : inline void ClosestPair2D::_set_label(Point * point, unsigned int review_flag) {
919 0 : if (point->review_flag == 0) _points_under_review.push_back(point);
920 0 : point->review_flag = review_flag;
921 0 : }
922 : void ClosestPair2D::remove(unsigned int ID) {
923 0 : Point * point_to_remove = & (_points[ID]);
924 0 : _remove_from_search_tree(point_to_remove);
925 0 : _deal_with_points_to_review();
926 0 : }
927 : void ClosestPair2D::_remove_from_search_tree(Point * point_to_remove) {
928 0 : _available_points.push(point_to_remove);
929 0 : _set_label(point_to_remove, _remove_heap_entry);
930 0 : unsigned int CP_range = min(_cp_search_range, size()-1);
931 0 : for (unsigned int ishift = 0; ishift < _nshift; ishift++) {
932 0 : circulator removed_circ = point_to_remove->circ[ishift];
933 0 : circulator right_end = removed_circ.next();
934 0 : _trees[ishift]->remove(removed_circ);
935 0 : circulator left_end = right_end, orig_right_end = right_end;
936 0 : for (unsigned int i = 0; i < CP_range; i++) {left_end--;}
937 0 : if (size()-1 < _cp_search_range) {
938 0 : left_end--; right_end--;
939 0 : }
940 : do {
941 0 : Point * left_point = left_end->point;
942 0 : if (left_point->neighbour == point_to_remove) {
943 : // we'll deal with it later...
944 0 : _add_label(left_point, _review_neighbour);
945 0 : } else {
946 : // check to see if right point has become its closest neighbour
947 0 : double dist2 = left_point->distance2(*right_end->point);
948 0 : if (dist2 < left_point->neighbour_dist2) {
949 0 : left_point->neighbour = right_end->point;
950 0 : left_point->neighbour_dist2 = dist2;
951 : // NB: (LESSER) REVIEW NEEDED HERE TOO...
952 0 : _add_label(left_point, _review_heap_entry);
953 0 : }
954 : }
955 0 : ++right_end;
956 0 : } while (++left_end != orig_right_end);
957 0 : } // ishift...
958 0 : }
959 : void ClosestPair2D::_deal_with_points_to_review() {
960 0 : unsigned int CP_range = min(_cp_search_range, size()-1);
961 0 : while(_points_under_review.size() > 0) {
962 0 : Point * this_point = _points_under_review.back();
963 0 : _points_under_review.pop_back();
964 0 : if (this_point->review_flag & _remove_heap_entry) {
965 0 : assert(!(this_point->review_flag ^ _remove_heap_entry));
966 0 : _heap->remove(_ID(this_point));
967 0 : }
968 : else {
969 0 : if (this_point->review_flag & _review_neighbour) {
970 0 : this_point->neighbour_dist2 = numeric_limits<double>::max();
971 : // among all three shifts
972 0 : for (unsigned int ishift = 0; ishift < _nshift; ishift++) {
973 0 : circulator other = this_point->circ[ishift];
974 : // among points within CP_range
975 0 : for (unsigned i=0; i < CP_range; i++) {
976 0 : ++other;
977 0 : double dist2 = this_point->distance2(*other->point);
978 0 : if (dist2 < this_point->neighbour_dist2) {
979 0 : this_point->neighbour_dist2 = dist2;
980 0 : this_point->neighbour = other->point;
981 0 : }
982 : }
983 0 : }
984 0 : }
985 0 : _heap->update(_ID(this_point), this_point->neighbour_dist2);
986 : }
987 0 : this_point->review_flag = 0;
988 : }
989 0 : }
990 : unsigned int ClosestPair2D::insert(const Coord2D & new_coord) {
991 0 : assert(_available_points.size() > 0);
992 0 : Point * new_point = _available_points.top();
993 0 : _available_points.pop();
994 0 : new_point->coord = new_coord;
995 0 : _insert_into_search_tree(new_point);
996 0 : _deal_with_points_to_review();
997 0 : return _ID(new_point);
998 : }
999 : unsigned int ClosestPair2D::replace(unsigned int ID1, unsigned int ID2,
1000 : const Coord2D & position) {
1001 0 : Point * point_to_remove = & (_points[ID1]);
1002 0 : _remove_from_search_tree(point_to_remove);
1003 0 : point_to_remove = & (_points[ID2]);
1004 0 : _remove_from_search_tree(point_to_remove);
1005 0 : Point * new_point = _available_points.top();
1006 0 : _available_points.pop();
1007 0 : new_point->coord = position;
1008 0 : _insert_into_search_tree(new_point);
1009 0 : _deal_with_points_to_review();
1010 0 : return _ID(new_point);
1011 : }
1012 : void ClosestPair2D::replace_many(
1013 : const std::vector<unsigned int> & IDs_to_remove,
1014 : const std::vector<Coord2D> & new_positions,
1015 : std::vector<unsigned int> & new_IDs) {
1016 0 : for (unsigned int i = 0; i < IDs_to_remove.size(); i++) {
1017 0 : _remove_from_search_tree(& (_points[IDs_to_remove[i]]));
1018 : }
1019 0 : new_IDs.resize(0);
1020 0 : for (unsigned int i = 0; i < new_positions.size(); i++) {
1021 0 : Point * new_point = _available_points.top();
1022 0 : _available_points.pop();
1023 0 : new_point->coord = new_positions[i];
1024 0 : _insert_into_search_tree(new_point);
1025 0 : new_IDs.push_back(_ID(new_point));
1026 : }
1027 0 : _deal_with_points_to_review();
1028 0 : }
1029 : void ClosestPair2D::_insert_into_search_tree(Point * new_point) {
1030 0 : _set_label(new_point, _review_heap_entry);
1031 0 : new_point->neighbour_dist2 = numeric_limits<double>::max();
1032 0 : unsigned int CP_range = min(_cp_search_range, size()-1);
1033 0 : for (unsigned ishift = 0; ishift < _nshift; ishift++) {
1034 0 : Shuffle new_shuffle;
1035 0 : _point2shuffle(*new_point, new_shuffle, _shifts[ishift]);
1036 0 : circulator new_circ = _trees[ishift]->insert(new_shuffle);
1037 0 : new_point->circ[ishift] = new_circ;
1038 0 : circulator right_edge = new_circ; right_edge++;
1039 0 : circulator left_edge = new_circ;
1040 0 : for (unsigned int i = 0; i < CP_range; i++) {left_edge--;}
1041 0 : do {
1042 0 : Point * left_point = left_edge->point;
1043 0 : Point * right_point = right_edge->point;
1044 0 : double new_dist2 = left_point->distance2(*new_point);
1045 0 : if (new_dist2 < left_point->neighbour_dist2) {
1046 0 : left_point->neighbour_dist2 = new_dist2;
1047 0 : left_point->neighbour = new_point;
1048 0 : _add_label(left_point, _review_heap_entry);
1049 0 : }
1050 0 : new_dist2 = new_point->distance2(*right_point);
1051 0 : if (new_dist2 < new_point->neighbour_dist2) {
1052 0 : new_point->neighbour_dist2 = new_dist2;
1053 0 : new_point->neighbour = right_point;
1054 0 : }
1055 0 : if (left_point->neighbour == right_point) {
1056 0 : _add_label(left_point, _review_neighbour);
1057 0 : }
1058 0 : right_edge++;
1059 0 : } while (++left_edge != new_circ);
1060 0 : }
1061 0 : }
1062 : FJCORE_END_NAMESPACE
1063 : #include<iostream>
1064 : #include<sstream>
1065 : #include<fstream>
1066 : #include<cmath>
1067 : #include<cstdlib>
1068 : #include<cassert>
1069 : #include<string>
1070 : #include<set>
1071 : FJCORE_BEGIN_NAMESPACE // defined in fastjet/internal/base.hh
1072 : using namespace std;
1073 : std::ostream * ClusterSequence::_fastjet_banner_ostr = &cout;
1074 0 : ClusterSequence::~ClusterSequence () {
1075 0 : if (_structure_shared_ptr()){
1076 0 : ClusterSequenceStructure* csi = dynamic_cast<ClusterSequenceStructure*>(_structure_shared_ptr());
1077 0 : assert(csi != NULL);
1078 0 : csi->set_associated_cs(NULL);
1079 0 : if (_deletes_self_when_unused) {
1080 0 : _structure_shared_ptr.set_count(_structure_shared_ptr.use_count()
1081 0 : + _structure_use_count_after_construction);
1082 : }
1083 0 : }
1084 0 : }
1085 : void ClusterSequence::signal_imminent_self_deletion() const {
1086 0 : assert(_deletes_self_when_unused);
1087 0 : _deletes_self_when_unused = false;
1088 0 : }
1089 : void ClusterSequence::_initialise_and_run (
1090 : const JetDefinition & jet_def_in,
1091 : const bool & writeout_combinations) {
1092 0 : _decant_options(jet_def_in, writeout_combinations);
1093 0 : _initialise_and_run_no_decant();
1094 0 : }
1095 : void ClusterSequence::_initialise_and_run_no_decant () {
1096 0 : _fill_initial_history();
1097 0 : if (n_particles() == 0) return;
1098 0 : if (_jet_algorithm == plugin_algorithm) {
1099 0 : _plugin_activated = true;
1100 0 : _jet_def.plugin()->run_clustering( (*this) );
1101 0 : _plugin_activated = false;
1102 0 : _update_structure_use_count();
1103 0 : return;
1104 0 : } else if (_jet_algorithm == ee_kt_algorithm ||
1105 0 : _jet_algorithm == ee_genkt_algorithm) {
1106 0 : _strategy = N2Plain;
1107 0 : if (_jet_algorithm == ee_kt_algorithm) {
1108 0 : assert(_Rparam > 2.0);
1109 0 : _invR2 = 1.0;
1110 0 : } else {
1111 0 : if (_Rparam > pi) {
1112 : // choose a value that ensures that back-to-back particles will
1113 : // always recombine
1114 : //_R2 = 4.0000000000001;
1115 0 : _R2 = 2 * ( 3.0 + cos(_Rparam) );
1116 0 : } else {
1117 0 : _R2 = 2 * ( 1.0 - cos(_Rparam) );
1118 : }
1119 0 : _invR2 = 1.0/_R2;
1120 : }
1121 0 : _simple_N2_cluster_EEBriefJet();
1122 0 : return;
1123 0 : } else if (_jet_algorithm == undefined_jet_algorithm) {
1124 0 : throw Error("A ClusterSequence cannot be created with an uninitialised JetDefinition");
1125 : }
1126 0 : if (_strategy == Best) {
1127 0 : int N = _jets.size();
1128 0 : if (min(1.0,max(0.1,_Rparam)*3.3)*N <= 30) {
1129 0 : _strategy = N2Plain;
1130 0 : } else if (N > 6200/pow(_Rparam,2.0) && _jet_def.jet_algorithm() == cambridge_algorithm) {
1131 0 : _strategy = NlnNCam;
1132 : #ifndef __FJCORE_DROP_CGAL
1133 : } else if ((N > 16000/pow(_Rparam,1.15) && _jet_def.jet_algorithm() != antikt_algorithm)
1134 : || N > 35000/pow(_Rparam,1.15)) {
1135 : _strategy = NlnN;
1136 : #endif // __FJCORE_DROP_CGAL
1137 0 : } else if (N <= 450) {
1138 0 : _strategy = N2Tiled;
1139 0 : } else {
1140 0 : _strategy = N2MinHeapTiled;
1141 : }
1142 0 : }
1143 0 : if (_Rparam >= twopi) {
1144 0 : if ( _strategy == NlnN
1145 0 : || _strategy == NlnN3pi
1146 0 : || _strategy == NlnNCam
1147 0 : || _strategy == NlnNCam2pi2R
1148 0 : || _strategy == NlnNCam4pi) {
1149 : #ifdef __FJCORE_DROP_CGAL
1150 0 : _strategy = N2MinHeapTiled;
1151 : #else
1152 : _strategy = NlnN4pi;
1153 : #endif
1154 0 : }
1155 0 : if (_jet_def.strategy() != Best && _strategy != _jet_def.strategy()) {
1156 0 : ostringstream oss;
1157 0 : oss << "Cluster strategy " << strategy_string(_jet_def.strategy())
1158 0 : << " automatically changed to " << strategy_string()
1159 0 : << " because the former is not supported for R = " << _Rparam
1160 0 : << " >= 2pi";
1161 0 : _changed_strategy_warning.warn(oss.str());
1162 0 : }
1163 : }
1164 0 : if (_strategy == N2Plain) {
1165 0 : this->_simple_N2_cluster_BriefJet();
1166 0 : } else if (_strategy == N2Tiled) {
1167 0 : this->_faster_tiled_N2_cluster();
1168 0 : } else if (_strategy == N2MinHeapTiled) {
1169 0 : this->_minheap_faster_tiled_N2_cluster();
1170 0 : } else if (_strategy == NlnN) {
1171 0 : this->_delaunay_cluster();
1172 0 : } else if (_strategy == NlnNCam) {
1173 0 : this->_CP2DChan_cluster_2piMultD();
1174 0 : } else if (_strategy == NlnN3pi || _strategy == NlnN4pi ) {
1175 0 : this->_delaunay_cluster();
1176 0 : } else if (_strategy == N3Dumb ) {
1177 0 : this->_really_dumb_cluster();
1178 0 : } else if (_strategy == N2PoorTiled) {
1179 0 : this->_tiled_N2_cluster();
1180 0 : } else if (_strategy == NlnNCam4pi) {
1181 0 : this->_CP2DChan_cluster();
1182 0 : } else if (_strategy == NlnNCam2pi2R) {
1183 0 : this->_CP2DChan_cluster_2pi2R();
1184 : } else {
1185 0 : ostringstream err;
1186 0 : err << "Unrecognised value for strategy: "<<_strategy;
1187 0 : throw Error(err.str());
1188 0 : }
1189 0 : }
1190 : bool ClusterSequence::_first_time = true;
1191 : int ClusterSequence::_n_exclusive_warnings = 0;
1192 : string fastjet_version_string() {
1193 0 : return "FastJet version "+string(fastjet_version)+" [fjcore]";
1194 0 : }
1195 : void ClusterSequence::print_banner() {
1196 0 : if (!_first_time) {return;}
1197 0 : _first_time = false;
1198 0 : ostream * ostr = _fastjet_banner_ostr;
1199 0 : if (!ostr) return;
1200 0 : (*ostr) << "#--------------------------------------------------------------------------\n";
1201 0 : (*ostr) << "# FastJet release " << fastjet_version << " [fjcore]" << endl;
1202 0 : (*ostr) << "# M. Cacciari, G.P. Salam and G. Soyez \n";
1203 0 : (*ostr) << "# A software package for jet finding and analysis at colliders \n";
1204 0 : (*ostr) << "# http://fastjet.fr \n";
1205 0 : (*ostr) << "# \n";
1206 0 : (*ostr) << "# Please cite EPJC72(2012)1896 [arXiv:1111.6097] if you use this package\n";
1207 0 : (*ostr) << "# for scientific work and optionally PLB641(2006)57 [hep-ph/0512210]. \n";
1208 0 : (*ostr) << "# \n";
1209 0 : (*ostr) << "# FastJet is provided without warranty under the terms of the GNU GPLv2.\n";
1210 0 : (*ostr) << "# It uses T. Chan's closest pair algorithm, S. Fortune's Voronoi code";
1211 : #ifndef __FJCORE_DROP_CGAL
1212 : (*ostr) << ",\n# CGAL ";
1213 : #else
1214 0 : (*ostr) << "\n# ";
1215 : #endif // __FJCORE_DROP_CGAL
1216 0 : (*ostr) << "and 3rd party plugin jet algorithms. See COPYING file for details.\n";
1217 0 : (*ostr) << "#--------------------------------------------------------------------------\n";
1218 0 : ostr->flush();
1219 0 : }
1220 : void ClusterSequence::_decant_options(const JetDefinition & jet_def_in,
1221 : const bool & writeout_combinations) {
1222 0 : _jet_def = jet_def_in;
1223 0 : _writeout_combinations = writeout_combinations;
1224 0 : _structure_shared_ptr.reset(new ClusterSequenceStructure(this));
1225 0 : _decant_options_partial();
1226 0 : }
1227 : void ClusterSequence::_decant_options_partial() {
1228 0 : print_banner();
1229 0 : _jet_algorithm = _jet_def.jet_algorithm();
1230 0 : _Rparam = _jet_def.R(); _R2 = _Rparam*_Rparam; _invR2 = 1.0/_R2;
1231 0 : _strategy = _jet_def.strategy();
1232 0 : _plugin_activated = false;
1233 0 : _update_structure_use_count(); // make sure it's correct already here
1234 0 : }
1235 : void ClusterSequence::_fill_initial_history () {
1236 0 : _jets.reserve(_jets.size()*2);
1237 0 : _history.reserve(_jets.size()*2);
1238 0 : _Qtot = 0;
1239 0 : for (int i = 0; i < static_cast<int>(_jets.size()) ; i++) {
1240 0 : history_element element;
1241 0 : element.parent1 = InexistentParent;
1242 0 : element.parent2 = InexistentParent;
1243 0 : element.child = Invalid;
1244 0 : element.jetp_index = i;
1245 0 : element.dij = 0.0;
1246 0 : element.max_dij_so_far = 0.0;
1247 0 : _history.push_back(element);
1248 0 : _jet_def.recombiner()->preprocess(_jets[i]);
1249 0 : _jets[i].set_cluster_hist_index(i);
1250 0 : _set_structure_shared_ptr(_jets[i]);
1251 0 : _Qtot += _jets[i].E();
1252 0 : }
1253 0 : _initial_n = _jets.size();
1254 0 : _deletes_self_when_unused = false;
1255 0 : }
1256 : string ClusterSequence::strategy_string (Strategy strategy_in) const {
1257 0 : string strategy;
1258 0 : switch(strategy_in) {
1259 : case NlnN:
1260 0 : strategy = "NlnN"; break;
1261 : case NlnN3pi:
1262 0 : strategy = "NlnN3pi"; break;
1263 : case NlnN4pi:
1264 0 : strategy = "NlnN4pi"; break;
1265 : case N2Plain:
1266 0 : strategy = "N2Plain"; break;
1267 : case N2Tiled:
1268 0 : strategy = "N2Tiled"; break;
1269 : case N2MinHeapTiled:
1270 0 : strategy = "N2MinHeapTiled"; break;
1271 : case N2PoorTiled:
1272 0 : strategy = "N2PoorTiled"; break;
1273 : case N3Dumb:
1274 0 : strategy = "N3Dumb"; break;
1275 : case NlnNCam4pi:
1276 0 : strategy = "NlnNCam4pi"; break;
1277 : case NlnNCam2pi2R:
1278 0 : strategy = "NlnNCam2pi2R"; break;
1279 : case NlnNCam:
1280 0 : strategy = "NlnNCam"; break; // 2piMultD
1281 : case plugin_strategy:
1282 0 : strategy = "plugin strategy"; break;
1283 : default:
1284 0 : strategy = "Unrecognized";
1285 : }
1286 : return strategy;
1287 0 : }
1288 : double ClusterSequence::jet_scale_for_algorithm(
1289 : const PseudoJet & jet) const {
1290 0 : if (_jet_algorithm == kt_algorithm) {return jet.kt2();}
1291 0 : else if (_jet_algorithm == cambridge_algorithm) {return 1.0;}
1292 0 : else if (_jet_algorithm == antikt_algorithm) {
1293 0 : double kt2=jet.kt2();
1294 0 : return kt2 > 1e-300 ? 1.0/kt2 : 1e300;
1295 0 : } else if (_jet_algorithm == genkt_algorithm) {
1296 0 : double kt2 = jet.kt2();
1297 0 : double p = jet_def().extra_param();
1298 0 : if (p <= 0 && kt2 < 1e-300) kt2 = 1e-300; // dodgy safety check
1299 0 : return pow(kt2, p);
1300 0 : } else if (_jet_algorithm == cambridge_for_passive_algorithm) {
1301 0 : double kt2 = jet.kt2();
1302 0 : double lim = _jet_def.extra_param();
1303 0 : if (kt2 < lim*lim && kt2 != 0.0) {
1304 0 : return 1.0/kt2;
1305 0 : } else {return 1.0;}
1306 0 : } else {throw Error("Unrecognised jet algorithm");}
1307 0 : }
1308 : void ClusterSequence::transfer_from_sequence(const ClusterSequence & from_seq,
1309 : const FunctionOfPseudoJet<PseudoJet> * action_on_jets){
1310 0 : if (will_delete_self_when_unused())
1311 0 : throw(Error("cannot use CS::transfer_from_sequence after a call to delete_self_when_unused()"));
1312 0 : _jet_def = from_seq._jet_def ;
1313 0 : _writeout_combinations = from_seq._writeout_combinations ;
1314 0 : _initial_n = from_seq._initial_n ;
1315 0 : _Rparam = from_seq._Rparam ;
1316 0 : _R2 = from_seq._R2 ;
1317 0 : _invR2 = from_seq._invR2 ;
1318 0 : _strategy = from_seq._strategy ;
1319 0 : _jet_algorithm = from_seq._jet_algorithm ;
1320 0 : _plugin_activated = from_seq._plugin_activated ;
1321 0 : if (action_on_jets)
1322 0 : _jets = (*action_on_jets)(from_seq._jets);
1323 : else
1324 0 : _jets = from_seq._jets;
1325 0 : _history = from_seq._history;
1326 0 : _extras = from_seq._extras;
1327 0 : if (_structure_shared_ptr()) {
1328 0 : if (_deletes_self_when_unused) throw Error("transfer_from_sequence cannot be used for a cluster sequence that deletes self when unused");
1329 0 : ClusterSequenceStructure* csi = dynamic_cast<ClusterSequenceStructure*>(_structure_shared_ptr());
1330 0 : assert(csi != NULL);
1331 0 : csi->set_associated_cs(NULL);
1332 0 : }
1333 0 : _structure_shared_ptr.reset(new ClusterSequenceStructure(this));
1334 0 : _update_structure_use_count();
1335 0 : for (unsigned int i=0; i<_jets.size(); i++){
1336 0 : _jets[i].set_cluster_hist_index(from_seq._jets[i].cluster_hist_index());
1337 0 : _set_structure_shared_ptr(_jets[i]);
1338 : }
1339 0 : }
1340 : void ClusterSequence::plugin_record_ij_recombination(
1341 : int jet_i, int jet_j, double dij,
1342 : const PseudoJet & newjet, int & newjet_k) {
1343 0 : plugin_record_ij_recombination(jet_i, jet_j, dij, newjet_k);
1344 0 : int tmp_index = _jets[newjet_k].cluster_hist_index();
1345 0 : _jets[newjet_k] = newjet;
1346 0 : _jets[newjet_k].set_cluster_hist_index(tmp_index);
1347 0 : _set_structure_shared_ptr(_jets[newjet_k]);
1348 0 : }
1349 : vector<PseudoJet> ClusterSequence::inclusive_jets (const double & ptmin) const{
1350 0 : double dcut = ptmin*ptmin;
1351 0 : int i = _history.size() - 1; // last jet
1352 0 : vector<PseudoJet> jets_local;
1353 0 : if (_jet_algorithm == kt_algorithm) {
1354 0 : while (i >= 0) {
1355 0 : if (_history[i].max_dij_so_far < dcut) {break;}
1356 0 : if (_history[i].parent2 == BeamJet && _history[i].dij >= dcut) {
1357 : // for beam jets
1358 0 : int parent1 = _history[i].parent1;
1359 0 : jets_local.push_back(_jets[_history[parent1].jetp_index]);}
1360 0 : i--;
1361 : }
1362 0 : } else if (_jet_algorithm == cambridge_algorithm) {
1363 0 : while (i >= 0) {
1364 0 : if (_history[i].parent2 != BeamJet) {break;}
1365 0 : int parent1 = _history[i].parent1;
1366 0 : const PseudoJet & jet = _jets[_history[parent1].jetp_index];
1367 0 : if (jet.perp2() >= dcut) {jets_local.push_back(jet);}
1368 0 : i--;
1369 : }
1370 0 : } else if (_jet_algorithm == plugin_algorithm
1371 0 : || _jet_algorithm == ee_kt_algorithm
1372 0 : || _jet_algorithm == antikt_algorithm
1373 0 : || _jet_algorithm == genkt_algorithm
1374 0 : || _jet_algorithm == ee_genkt_algorithm
1375 0 : || _jet_algorithm == cambridge_for_passive_algorithm) {
1376 0 : while (i >= 0) {
1377 0 : if (_history[i].parent2 == BeamJet) {
1378 0 : int parent1 = _history[i].parent1;
1379 0 : const PseudoJet & jet = _jets[_history[parent1].jetp_index];
1380 0 : if (jet.perp2() >= dcut) {jets_local.push_back(jet);}
1381 0 : }
1382 0 : i--;
1383 : }
1384 0 : } else {throw Error("cs::inclusive_jets(...): Unrecognized jet algorithm");}
1385 : return jets_local;
1386 0 : }
1387 : int ClusterSequence::n_exclusive_jets (const double & dcut) const {
1388 0 : int i = _history.size() - 1; // last jet
1389 0 : while (i >= 0) {
1390 0 : if (_history[i].max_dij_so_far <= dcut) {break;}
1391 0 : i--;
1392 : }
1393 0 : int stop_point = i + 1;
1394 0 : int njets = 2*_initial_n - stop_point;
1395 0 : return njets;
1396 : }
1397 : vector<PseudoJet> ClusterSequence::exclusive_jets (const double & dcut) const {
1398 0 : int njets = n_exclusive_jets(dcut);
1399 0 : return exclusive_jets(njets);
1400 0 : }
1401 : vector<PseudoJet> ClusterSequence::exclusive_jets (const int & njets) const {
1402 0 : if (njets > _initial_n) {
1403 0 : ostringstream err;
1404 0 : err << "Requested " << njets << " exclusive jets, but there were only "
1405 0 : << _initial_n << " particles in the event";
1406 0 : throw Error(err.str());
1407 0 : }
1408 0 : return exclusive_jets_up_to(njets);
1409 0 : }
1410 : vector<PseudoJet> ClusterSequence::exclusive_jets_up_to (const int & njets) const {
1411 0 : if (( _jet_def.jet_algorithm() != kt_algorithm) &&
1412 0 : ( _jet_def.jet_algorithm() != cambridge_algorithm) &&
1413 0 : ( _jet_def.jet_algorithm() != ee_kt_algorithm) &&
1414 0 : (((_jet_def.jet_algorithm() != genkt_algorithm) &&
1415 0 : (_jet_def.jet_algorithm() != ee_genkt_algorithm)) ||
1416 0 : (_jet_def.extra_param() <0)) &&
1417 0 : ((_jet_def.jet_algorithm() != plugin_algorithm) ||
1418 0 : (!_jet_def.plugin()->exclusive_sequence_meaningful())) &&
1419 0 : (_n_exclusive_warnings < 5)) {
1420 0 : _n_exclusive_warnings++;
1421 0 : cerr << "FastJet WARNING: dcut and exclusive jets for jet-finders other than kt should be interpreted with care." << endl;
1422 0 : }
1423 0 : int stop_point = 2*_initial_n - njets;
1424 0 : if (stop_point < _initial_n) stop_point = _initial_n;
1425 0 : if (2*_initial_n != static_cast<int>(_history.size())) {
1426 0 : ostringstream err;
1427 0 : err << "2*_initial_n != _history.size() -- this endangers internal assumptions!\n";
1428 0 : throw Error(err.str());
1429 0 : }
1430 0 : vector<PseudoJet> jets_local;
1431 0 : for (unsigned int i = stop_point; i < _history.size(); i++) {
1432 0 : int parent1 = _history[i].parent1;
1433 0 : if (parent1 < stop_point) {
1434 0 : jets_local.push_back(_jets[_history[parent1].jetp_index]);
1435 : }
1436 0 : int parent2 = _history[i].parent2;
1437 0 : if (parent2 < stop_point && parent2 > 0) {
1438 0 : jets_local.push_back(_jets[_history[parent2].jetp_index]);
1439 : }
1440 : }
1441 0 : if (int(jets_local.size()) != min(_initial_n, njets)) {
1442 0 : ostringstream err;
1443 0 : err << "ClusterSequence::exclusive_jets: size of returned vector ("
1444 0 : <<jets_local.size()<<") does not coincide with requested number of jets ("
1445 0 : <<njets<<")";
1446 0 : throw Error(err.str());
1447 0 : }
1448 : return jets_local;
1449 0 : }
1450 : double ClusterSequence::exclusive_dmerge (const int & njets) const {
1451 0 : assert(njets >= 0);
1452 0 : if (njets >= _initial_n) {return 0.0;}
1453 0 : return _history[2*_initial_n-njets-1].dij;
1454 0 : }
1455 : double ClusterSequence::exclusive_dmerge_max (const int & njets) const {
1456 0 : assert(njets >= 0);
1457 0 : if (njets >= _initial_n) {return 0.0;}
1458 0 : return _history[2*_initial_n-njets-1].max_dij_so_far;
1459 0 : }
1460 : std::vector<PseudoJet> ClusterSequence::exclusive_subjets
1461 : (const PseudoJet & jet, const double & dcut) const {
1462 0 : set<const history_element*> subhist;
1463 0 : get_subhist_set(subhist, jet, dcut, 0);
1464 0 : vector<PseudoJet> subjets;
1465 0 : subjets.reserve(subhist.size());
1466 0 : for (set<const history_element*>::iterator elem = subhist.begin();
1467 0 : elem != subhist.end(); elem++) {
1468 0 : subjets.push_back(_jets[(*elem)->jetp_index]);
1469 : }
1470 : return subjets;
1471 0 : }
1472 : int ClusterSequence::n_exclusive_subjets(const PseudoJet & jet,
1473 : const double & dcut) const {
1474 0 : set<const history_element*> subhist;
1475 0 : get_subhist_set(subhist, jet, dcut, 0);
1476 0 : return subhist.size();
1477 0 : }
1478 : std::vector<PseudoJet> ClusterSequence::exclusive_subjets
1479 : (const PseudoJet & jet, int nsub) const {
1480 0 : vector<PseudoJet> subjets = exclusive_subjets_up_to(jet, nsub);
1481 0 : if (int(subjets.size()) < nsub) {
1482 0 : ostringstream err;
1483 0 : err << "Requested " << nsub << " exclusive subjets, but there were only "
1484 0 : << subjets.size() << " particles in the jet";
1485 0 : throw Error(err.str());
1486 0 : }
1487 : return subjets;
1488 0 : }
1489 : std::vector<PseudoJet> ClusterSequence::exclusive_subjets_up_to
1490 : (const PseudoJet & jet, int nsub) const {
1491 0 : set<const history_element*> subhist;
1492 0 : vector<PseudoJet> subjets;
1493 0 : if (nsub < 0) throw Error("Requested a negative number of subjets. This is nonsensical.");
1494 0 : if (nsub == 0) return subjets;
1495 0 : get_subhist_set(subhist, jet, -1.0, nsub);
1496 0 : subjets.reserve(subhist.size());
1497 0 : for (set<const history_element*>::iterator elem = subhist.begin();
1498 0 : elem != subhist.end(); elem++) {
1499 0 : subjets.push_back(_jets[(*elem)->jetp_index]);
1500 : }
1501 0 : return subjets;
1502 0 : }
1503 : double ClusterSequence::exclusive_subdmerge(const PseudoJet & jet, int nsub) const {
1504 0 : set<const history_element*> subhist;
1505 0 : get_subhist_set(subhist, jet, -1.0, nsub);
1506 0 : set<const history_element*>::iterator highest = subhist.end();
1507 0 : highest--;
1508 0 : return (*highest)->dij;
1509 0 : }
1510 : double ClusterSequence::exclusive_subdmerge_max(const PseudoJet & jet, int nsub) const {
1511 0 : set<const history_element*> subhist;
1512 0 : get_subhist_set(subhist, jet, -1.0, nsub);
1513 0 : set<const history_element*>::iterator highest = subhist.end();
1514 0 : highest--;
1515 0 : return (*highest)->max_dij_so_far;
1516 0 : }
1517 : void ClusterSequence::get_subhist_set(set<const history_element*> & subhist,
1518 : const PseudoJet & jet,
1519 : double dcut, int maxjet) const {
1520 0 : assert(contains(jet));
1521 0 : subhist.clear();
1522 0 : subhist.insert(&(_history[jet.cluster_hist_index()]));
1523 : int njet = 1;
1524 0 : while (true) {
1525 0 : set<const history_element*>::iterator highest = subhist.end();
1526 0 : assert (highest != subhist.begin());
1527 0 : highest--;
1528 0 : const history_element* elem = *highest;
1529 0 : if (njet == maxjet) break;
1530 0 : if (elem->parent1 < 0) break;
1531 0 : if (elem->max_dij_so_far <= dcut) break;
1532 0 : subhist.erase(highest);
1533 0 : subhist.insert(&(_history[elem->parent1]));
1534 0 : subhist.insert(&(_history[elem->parent2]));
1535 0 : njet++;
1536 0 : }
1537 0 : }
1538 : bool ClusterSequence::object_in_jet(const PseudoJet & object,
1539 : const PseudoJet & jet) const {
1540 0 : assert(contains(object) && contains(jet));
1541 : const PseudoJet * this_object = &object;
1542 0 : const PseudoJet * childp;
1543 0 : while(true) {
1544 0 : if (this_object->cluster_hist_index() == jet.cluster_hist_index()) {
1545 0 : return true;
1546 0 : } else if (has_child(*this_object, childp)) {
1547 0 : this_object = childp;
1548 : } else {
1549 0 : return false;
1550 : }
1551 : }
1552 0 : }
1553 : bool ClusterSequence::has_parents(const PseudoJet & jet, PseudoJet & parent1,
1554 : PseudoJet & parent2) const {
1555 0 : const history_element & hist = _history[jet.cluster_hist_index()];
1556 0 : assert ((hist.parent1 >= 0 && hist.parent2 >= 0) ||
1557 : (hist.parent1 < 0 && hist.parent2 < 0));
1558 0 : if (hist.parent1 < 0) {
1559 0 : parent1 = PseudoJet(0.0,0.0,0.0,0.0);
1560 0 : parent2 = parent1;
1561 0 : return false;
1562 : } else {
1563 0 : parent1 = _jets[_history[hist.parent1].jetp_index];
1564 0 : parent2 = _jets[_history[hist.parent2].jetp_index];
1565 0 : if (parent1.perp2() < parent2.perp2()) std::swap(parent1,parent2);
1566 0 : return true;
1567 : }
1568 0 : }
1569 : bool ClusterSequence::has_child(const PseudoJet & jet, PseudoJet & child) const {
1570 0 : const PseudoJet * childp;
1571 0 : bool res = has_child(jet, childp);
1572 0 : if (res) {
1573 0 : child = *childp;
1574 0 : return true;
1575 : } else {
1576 0 : child = PseudoJet(0.0,0.0,0.0,0.0);
1577 0 : return false;
1578 : }
1579 0 : }
1580 : bool ClusterSequence::has_child(const PseudoJet & jet, const PseudoJet * & childp) const {
1581 0 : const history_element & hist = _history[jet.cluster_hist_index()];
1582 0 : if (hist.child >= 0 && _history[hist.child].jetp_index >= 0) {
1583 0 : childp = &(_jets[_history[hist.child].jetp_index]);
1584 0 : return true;
1585 : } else {
1586 0 : childp = NULL;
1587 0 : return false;
1588 : }
1589 0 : }
1590 : bool ClusterSequence::has_partner(const PseudoJet & jet,
1591 : PseudoJet & partner) const {
1592 0 : const history_element & hist = _history[jet.cluster_hist_index()];
1593 0 : if (hist.child >= 0 && _history[hist.child].parent2 >= 0) {
1594 0 : const history_element & child_hist = _history[hist.child];
1595 0 : if (child_hist.parent1 == jet.cluster_hist_index()) {
1596 0 : partner = _jets[_history[child_hist.parent2].jetp_index];
1597 0 : } else {
1598 0 : partner = _jets[_history[child_hist.parent1].jetp_index];
1599 : }
1600 : return true;
1601 : } else {
1602 0 : partner = PseudoJet(0.0,0.0,0.0,0.0);
1603 0 : return false;
1604 : }
1605 0 : }
1606 : vector<PseudoJet> ClusterSequence::constituents (const PseudoJet & jet) const {
1607 0 : vector<PseudoJet> subjets;
1608 0 : add_constituents(jet, subjets);
1609 : return subjets;
1610 0 : }
1611 : void ClusterSequence::print_jets_for_root(const std::vector<PseudoJet> & jets_in,
1612 : ostream & ostr) const {
1613 0 : for (unsigned i = 0; i < jets_in.size(); i++) {
1614 0 : ostr << i << " "
1615 0 : << jets_in[i].px() << " "
1616 0 : << jets_in[i].py() << " "
1617 0 : << jets_in[i].pz() << " "
1618 0 : << jets_in[i].E() << endl;
1619 0 : vector<PseudoJet> cst = constituents(jets_in[i]);
1620 0 : for (unsigned j = 0; j < cst.size() ; j++) {
1621 0 : ostr << " " << j << " "
1622 0 : << cst[j].rap() << " "
1623 0 : << cst[j].phi() << " "
1624 0 : << cst[j].perp() << endl;
1625 : }
1626 0 : ostr << "#END" << endl;
1627 0 : }
1628 0 : }
1629 : void ClusterSequence::print_jets_for_root(const std::vector<PseudoJet> & jets_in,
1630 : const std::string & filename,
1631 : const std::string & comment ) const {
1632 0 : std::ofstream ostr(filename.c_str());
1633 0 : if (comment != "") ostr << "# " << comment << endl;
1634 0 : print_jets_for_root(jets_in, ostr);
1635 0 : }
1636 : vector<int> ClusterSequence::particle_jet_indices(
1637 : const vector<PseudoJet> & jets_in) const {
1638 0 : vector<int> indices(n_particles());
1639 0 : for (unsigned ipart = 0; ipart < n_particles(); ipart++)
1640 0 : indices[ipart] = -1;
1641 0 : for (unsigned ijet = 0; ijet < jets_in.size(); ijet++) {
1642 0 : vector<PseudoJet> jet_constituents(constituents(jets_in[ijet]));
1643 0 : for (unsigned ip = 0; ip < jet_constituents.size(); ip++) {
1644 0 : unsigned iclust = jet_constituents[ip].cluster_hist_index();
1645 0 : unsigned ipart = history()[iclust].jetp_index;
1646 0 : indices[ipart] = ijet;
1647 : }
1648 0 : }
1649 : return indices;
1650 0 : }
1651 : void ClusterSequence::add_constituents (
1652 : const PseudoJet & jet, vector<PseudoJet> & subjet_vector) const {
1653 0 : int i = jet.cluster_hist_index();
1654 0 : int parent1 = _history[i].parent1;
1655 0 : int parent2 = _history[i].parent2;
1656 0 : if (parent1 == InexistentParent) {
1657 0 : subjet_vector.push_back(_jets[i]);
1658 0 : return;
1659 : }
1660 0 : add_constituents(_jets[_history[parent1].jetp_index], subjet_vector);
1661 0 : if (parent2 != BeamJet) {
1662 0 : add_constituents(_jets[_history[parent2].jetp_index], subjet_vector);
1663 0 : }
1664 0 : }
1665 : void ClusterSequence::_add_step_to_history (
1666 : const int & step_number, const int & parent1,
1667 : const int & parent2, const int & jetp_index,
1668 : const double & dij) {
1669 0 : history_element element;
1670 0 : element.parent1 = parent1;
1671 0 : element.parent2 = parent2;
1672 0 : element.jetp_index = jetp_index;
1673 0 : element.child = Invalid;
1674 0 : element.dij = dij;
1675 0 : element.max_dij_so_far = max(dij,_history[_history.size()-1].max_dij_so_far);
1676 0 : _history.push_back(element);
1677 0 : int local_step = _history.size()-1;
1678 0 : assert(local_step == step_number);
1679 0 : assert(parent1 >= 0);
1680 0 : _history[parent1].child = local_step;
1681 0 : if (parent2 >= 0) {_history[parent2].child = local_step;}
1682 0 : if (jetp_index != Invalid) {
1683 0 : assert(jetp_index >= 0);
1684 0 : _jets[jetp_index].set_cluster_hist_index(local_step);
1685 0 : _set_structure_shared_ptr(_jets[jetp_index]);
1686 0 : }
1687 0 : if (_writeout_combinations) {
1688 0 : cout << local_step << ": "
1689 0 : << parent1 << " with " << parent2
1690 0 : << "; y = "<< dij<<endl;
1691 0 : }
1692 0 : }
1693 : vector<int> ClusterSequence::unique_history_order() const {
1694 0 : valarray<int> lowest_constituent(_history.size());
1695 0 : int hist_n = _history.size();
1696 0 : lowest_constituent = hist_n; // give it a large number
1697 0 : for (int i = 0; i < hist_n; i++) {
1698 0 : lowest_constituent[i] = min(lowest_constituent[i],i);
1699 0 : if (_history[i].child > 0) lowest_constituent[_history[i].child]
1700 0 : = min(lowest_constituent[_history[i].child],lowest_constituent[i]);
1701 : }
1702 0 : valarray<bool> extracted(_history.size()); extracted = false;
1703 0 : vector<int> unique_tree;
1704 0 : unique_tree.reserve(_history.size());
1705 0 : for (unsigned i = 0; i < n_particles(); i++) {
1706 0 : if (!extracted[i]) {
1707 0 : unique_tree.push_back(i);
1708 0 : extracted[i] = true;
1709 0 : _extract_tree_children(i, extracted, lowest_constituent, unique_tree);
1710 : }
1711 : }
1712 : return unique_tree;
1713 0 : }
1714 : void ClusterSequence::_extract_tree_children(
1715 : int position,
1716 : valarray<bool> & extracted,
1717 : const valarray<int> & lowest_constituent,
1718 : vector<int> & unique_tree) const {
1719 0 : if (!extracted[position]) {
1720 0 : _extract_tree_parents(position,extracted,lowest_constituent,unique_tree);
1721 0 : }
1722 0 : int child = _history[position].child;
1723 0 : if (child >= 0) _extract_tree_children(child,extracted,lowest_constituent,unique_tree);
1724 0 : }
1725 : vector<PseudoJet> ClusterSequence::unclustered_particles() const {
1726 0 : vector<PseudoJet> unclustered;
1727 0 : for (unsigned i = 0; i < n_particles() ; i++) {
1728 0 : if (_history[i].child == Invalid)
1729 0 : unclustered.push_back(_jets[_history[i].jetp_index]);
1730 : }
1731 : return unclustered;
1732 0 : }
1733 : vector<PseudoJet> ClusterSequence::childless_pseudojets() const {
1734 0 : vector<PseudoJet> unclustered;
1735 0 : for (unsigned i = 0; i < _history.size() ; i++) {
1736 0 : if ((_history[i].child == Invalid) && (_history[i].parent2 != BeamJet))
1737 0 : unclustered.push_back(_jets[_history[i].jetp_index]);
1738 : }
1739 : return unclustered;
1740 0 : }
1741 : bool ClusterSequence::contains(const PseudoJet & jet) const {
1742 0 : return jet.cluster_hist_index() >= 0
1743 0 : && jet.cluster_hist_index() < int(_history.size())
1744 0 : && jet.has_valid_cluster_sequence()
1745 0 : && jet.associated_cluster_sequence() == this;
1746 : }
1747 : void ClusterSequence::_extract_tree_parents(
1748 : int position,
1749 : valarray<bool> & extracted,
1750 : const valarray<int> & lowest_constituent,
1751 : vector<int> & unique_tree) const {
1752 0 : if (!extracted[position]) {
1753 0 : int parent1 = _history[position].parent1;
1754 0 : int parent2 = _history[position].parent2;
1755 0 : if (parent1 >= 0 && parent2 >= 0) {
1756 0 : if (lowest_constituent[parent1] > lowest_constituent[parent2])
1757 0 : std::swap(parent1, parent2);
1758 : }
1759 0 : if (parent1 >= 0 && !extracted[parent1])
1760 0 : _extract_tree_parents(parent1,extracted,lowest_constituent,unique_tree);
1761 0 : if (parent2 >= 0 && !extracted[parent2])
1762 0 : _extract_tree_parents(parent2,extracted,lowest_constituent,unique_tree);
1763 0 : unique_tree.push_back(position);
1764 0 : extracted[position] = true;
1765 0 : }
1766 0 : }
1767 : void ClusterSequence::_do_ij_recombination_step(
1768 : const int & jet_i, const int & jet_j,
1769 : const double & dij,
1770 : int & newjet_k) {
1771 0 : PseudoJet newjet(false);
1772 0 : _jet_def.recombiner()->recombine(_jets[jet_i], _jets[jet_j], newjet);
1773 0 : _jets.push_back(newjet);
1774 0 : newjet_k = _jets.size()-1;
1775 0 : int newstep_k = _history.size();
1776 0 : _jets[newjet_k].set_cluster_hist_index(newstep_k);
1777 0 : int hist_i = _jets[jet_i].cluster_hist_index();
1778 0 : int hist_j = _jets[jet_j].cluster_hist_index();
1779 0 : _add_step_to_history(newstep_k, min(hist_i, hist_j), max(hist_i,hist_j),
1780 : newjet_k, dij);
1781 0 : }
1782 : void ClusterSequence::_do_iB_recombination_step(
1783 : const int & jet_i, const double & diB) {
1784 0 : int newstep_k = _history.size();
1785 0 : _add_step_to_history(newstep_k,_jets[jet_i].cluster_hist_index(),BeamJet,
1786 0 : Invalid, diB);
1787 0 : }
1788 6 : LimitedWarning ClusterSequence::_changed_strategy_warning;
1789 : void ClusterSequence::_set_structure_shared_ptr(PseudoJet & j) {
1790 0 : j.set_structure_shared_ptr(_structure_shared_ptr);
1791 0 : _update_structure_use_count();
1792 0 : }
1793 : void ClusterSequence::_update_structure_use_count() {
1794 0 : _structure_use_count_after_construction = _structure_shared_ptr.use_count();
1795 0 : }
1796 : void ClusterSequence::delete_self_when_unused() {
1797 0 : int new_count = _structure_shared_ptr.use_count() - _structure_use_count_after_construction;
1798 0 : if (new_count <= 0) {
1799 0 : 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");
1800 : }
1801 0 : _structure_shared_ptr.set_count(new_count);
1802 0 : _deletes_self_when_unused = true;
1803 0 : }
1804 : FJCORE_END_NAMESPACE
1805 : #include<limits>
1806 : #include<vector>
1807 : #include<cmath>
1808 : #include<iostream>
1809 : FJCORE_BEGIN_NAMESPACE // defined in fastjet/internal/base.hh
1810 : using namespace std;
1811 : namespace Private {
1812 : class MirrorInfo{
1813 : public:
1814 : int orig, mirror;
1815 0 : MirrorInfo(int a, int b) : orig(a), mirror(b) {};
1816 0 : MirrorInfo() {};
1817 : };
1818 : bool make_mirror(Coord2D & point, double Dlim) {
1819 0 : if (point.y < Dlim) {point.y += twopi; return true;}
1820 0 : if (twopi-point.y < Dlim) {point.y -= twopi; return true;}
1821 0 : return false;
1822 0 : }
1823 : }
1824 : using namespace Private;
1825 : void ClusterSequence::_CP2DChan_limited_cluster (double Dlim) {
1826 0 : unsigned int n = _initial_n;
1827 0 : vector<MirrorInfo> coordIDs(2*n); // coord IDs of a given jetID
1828 0 : vector<int> jetIDs(2*n); // jet ID for a given coord ID
1829 0 : vector<Coord2D> coords(2*n); // our coordinates (and copies)
1830 0 : double Dlim4mirror = min(Dlim,pi);
1831 0 : double minrap = numeric_limits<double>::max();
1832 0 : double maxrap = -minrap;
1833 : int coord_index = -1;
1834 : int n_active = 0;
1835 0 : for (unsigned jet_i = 0; jet_i < _jets.size(); jet_i++) {
1836 0 : if (_history[_jets[jet_i].cluster_hist_index()].child != Invalid ||
1837 0 : (_jets[jet_i].E() == abs(_jets[jet_i].pz()) &&
1838 0 : _jets[jet_i].perp2() == 0.0)
1839 : ) {continue;}
1840 0 : n_active++;
1841 0 : coordIDs[jet_i].orig = ++coord_index;
1842 0 : coords[coord_index] = Coord2D(_jets[jet_i].rap(), _jets[jet_i].phi_02pi());
1843 0 : jetIDs[coord_index] = jet_i;
1844 0 : minrap = min(coords[coord_index].x,minrap);
1845 0 : maxrap = max(coords[coord_index].x,maxrap);
1846 0 : Coord2D mirror_point(coords[coord_index]);
1847 0 : if (make_mirror(mirror_point, Dlim4mirror)) {
1848 0 : coordIDs[jet_i].mirror = ++coord_index;
1849 0 : coords[coord_index] = mirror_point;
1850 0 : jetIDs[coord_index] = jet_i;
1851 0 : } else {
1852 0 : coordIDs[jet_i].mirror = Invalid;
1853 : }
1854 0 : }
1855 0 : coords.resize(coord_index+1);
1856 0 : Coord2D left_edge(minrap-1.0, -3.15); // a security margin below -pi
1857 0 : Coord2D right_edge(maxrap+1.0, 9.45); // a security margin above 3*pi
1858 0 : ClosestPair2D cp(coords, left_edge, right_edge);
1859 0 : vector<Coord2D> new_points(2);
1860 0 : vector<unsigned int> cIDs_to_remove(4);
1861 0 : vector<unsigned int> new_cIDs(2);
1862 : do {
1863 0 : unsigned int cID1, cID2;
1864 0 : double distance2;
1865 0 : cp.closest_pair(cID1,cID2,distance2);
1866 0 : if (distance2 > Dlim*Dlim) {break;}
1867 0 : distance2 *= _invR2;
1868 0 : int jet_i = jetIDs[cID1];
1869 0 : int jet_j = jetIDs[cID2];
1870 0 : assert (jet_i != jet_j); // to catch issue of recombining with mirror point
1871 0 : int newjet_k;
1872 0 : _do_ij_recombination_step(jet_i, jet_j, distance2, newjet_k);
1873 0 : if (--n_active == 1) {break;}
1874 0 : cIDs_to_remove.resize(0);
1875 0 : cIDs_to_remove.push_back(coordIDs[jet_i].orig);
1876 0 : cIDs_to_remove.push_back(coordIDs[jet_j].orig);
1877 0 : if (coordIDs[jet_i].mirror != Invalid)
1878 0 : cIDs_to_remove.push_back(coordIDs[jet_i].mirror);
1879 0 : if (coordIDs[jet_j].mirror != Invalid)
1880 0 : cIDs_to_remove.push_back(coordIDs[jet_j].mirror);
1881 0 : Coord2D new_point(_jets[newjet_k].rap(),_jets[newjet_k].phi_02pi());
1882 0 : new_points.resize(0);
1883 0 : new_points.push_back(new_point);
1884 0 : if (make_mirror(new_point, Dlim4mirror)) new_points.push_back(new_point); //< same warning as before concerning the mirroring
1885 0 : cp.replace_many(cIDs_to_remove, new_points, new_cIDs);
1886 0 : coordIDs[newjet_k].orig = new_cIDs[0];
1887 0 : jetIDs[new_cIDs[0]] = newjet_k;
1888 0 : if (new_cIDs.size() == 2) {
1889 0 : coordIDs[newjet_k].mirror = new_cIDs[1];
1890 0 : jetIDs[new_cIDs[1]] = newjet_k;
1891 0 : } else {coordIDs[newjet_k].mirror = Invalid;}
1892 0 : } while(true);
1893 0 : }
1894 : void ClusterSequence::_CP2DChan_cluster_2pi2R () {
1895 0 : if (_jet_algorithm != cambridge_algorithm) throw Error("CP2DChan clustering method called for a jet-finder that is not the cambridge algorithm");
1896 0 : _CP2DChan_limited_cluster(_Rparam);
1897 0 : _do_Cambridge_inclusive_jets();
1898 0 : }
1899 : void ClusterSequence::_CP2DChan_cluster_2piMultD () {
1900 0 : if (_Rparam >= 0.39) {
1901 0 : _CP2DChan_limited_cluster(min(_Rparam/2,0.3));
1902 0 : }
1903 0 : _CP2DChan_cluster_2pi2R ();
1904 0 : }
1905 : void ClusterSequence::_CP2DChan_cluster () {
1906 0 : if (_jet_algorithm != cambridge_algorithm) throw Error("_CP2DChan_cluster called for a jet-finder that is not the cambridge algorithm");
1907 0 : unsigned int n = _jets.size();
1908 0 : vector<MirrorInfo> coordIDs(2*n); // link from original to mirror indices
1909 0 : vector<int> jetIDs(2*n); // link from mirror to original indices
1910 0 : vector<Coord2D> coords(2*n); // our coordinates (and copies)
1911 0 : double minrap = numeric_limits<double>::max();
1912 0 : double maxrap = -minrap;
1913 : int coord_index = 0;
1914 0 : for (unsigned i = 0; i < n; i++) {
1915 0 : if (_jets[i].E() == abs(_jets[i].pz()) && _jets[i].perp2() == 0.0) {
1916 0 : coordIDs[i] = MirrorInfo(BeamJet,BeamJet);
1917 0 : } else {
1918 0 : coordIDs[i].orig = coord_index;
1919 0 : coordIDs[i].mirror = coord_index+1;
1920 0 : coords[coord_index] = Coord2D(_jets[i].rap(), _jets[i].phi_02pi());
1921 0 : coords[coord_index+1] = Coord2D(_jets[i].rap(), _jets[i].phi_02pi()+twopi);
1922 0 : jetIDs[coord_index] = i;
1923 0 : jetIDs[coord_index+1] = i;
1924 0 : minrap = min(coords[coord_index].x,minrap);
1925 0 : maxrap = max(coords[coord_index].x,maxrap);
1926 0 : coord_index += 2;
1927 : }
1928 : }
1929 0 : for (unsigned i = n; i < 2*n; i++) {coordIDs[i].orig = Invalid;}
1930 0 : coords.resize(coord_index);
1931 0 : Coord2D left_edge(minrap-1.0, 0.0);
1932 0 : Coord2D right_edge(maxrap+1.0, 2*twopi);
1933 0 : ClosestPair2D cp(coords, left_edge, right_edge);
1934 0 : vector<Coord2D> new_points(2);
1935 0 : vector<unsigned int> cIDs_to_remove(4);
1936 0 : vector<unsigned int> new_cIDs(2);
1937 : do {
1938 0 : unsigned int cID1, cID2;
1939 0 : double distance2;
1940 0 : cp.closest_pair(cID1,cID2,distance2);
1941 0 : distance2 *= _invR2;
1942 0 : if (distance2 > 1.0) {break;}
1943 0 : int jet_i = jetIDs[cID1];
1944 0 : int jet_j = jetIDs[cID2];
1945 0 : assert (jet_i != jet_j); // to catch issue of recombining with mirror point
1946 0 : int newjet_k;
1947 0 : _do_ij_recombination_step(jet_i, jet_j, distance2, newjet_k);
1948 0 : cIDs_to_remove[0] = coordIDs[jet_i].orig;
1949 0 : cIDs_to_remove[1] = coordIDs[jet_i].mirror;
1950 0 : cIDs_to_remove[2] = coordIDs[jet_j].orig;
1951 0 : cIDs_to_remove[3] = coordIDs[jet_j].mirror;
1952 0 : new_points[0] = Coord2D(_jets[newjet_k].rap(),_jets[newjet_k].phi_02pi());
1953 0 : new_points[1] = Coord2D(_jets[newjet_k].rap(),_jets[newjet_k].phi_02pi()+twopi);
1954 0 : new_cIDs[0] = cp.replace(cIDs_to_remove[0], cIDs_to_remove[2], new_points[0]);
1955 0 : new_cIDs[1] = cp.replace(cIDs_to_remove[1], cIDs_to_remove[3], new_points[1]);
1956 0 : coordIDs[jet_i].orig = Invalid;
1957 0 : coordIDs[jet_j].orig = Invalid;
1958 0 : coordIDs[newjet_k] = MirrorInfo(new_cIDs[0], new_cIDs[1]);
1959 0 : jetIDs[new_cIDs[0]] = newjet_k;
1960 0 : jetIDs[new_cIDs[1]] = newjet_k;
1961 0 : n--;
1962 0 : if (n == 1) {break;}
1963 0 : } while(true);
1964 0 : _do_Cambridge_inclusive_jets();
1965 0 : }
1966 : void ClusterSequence::_do_Cambridge_inclusive_jets () {
1967 0 : unsigned int n = _history.size();
1968 0 : for (unsigned int hist_i = 0; hist_i < n; hist_i++) {
1969 0 : if (_history[hist_i].child == Invalid) {
1970 0 : _do_iB_recombination_step(_history[hist_i].jetp_index, 1.0);
1971 0 : }
1972 : }
1973 0 : }
1974 : FJCORE_END_NAMESPACE
1975 : #include<iostream>
1976 : #include<sstream>
1977 : #include<cmath>
1978 : #include <cstdlib>
1979 : #include<cassert>
1980 : #include<memory>
1981 : #ifndef __FJCORE_DROP_CGAL // in case we do not have the code for CGAL
1982 : #include "fastjet/internal/Dnn4piCylinder.hh"
1983 : #include "fastjet/internal/Dnn3piCylinder.hh"
1984 : #include "fastjet/internal/Dnn2piCylinder.hh"
1985 : #endif // __FJCORE_DROP_CGAL
1986 : FJCORE_BEGIN_NAMESPACE // defined in fastjet/internal/base.hh
1987 : using namespace std;
1988 : void ClusterSequence::_delaunay_cluster () {
1989 0 : int n = _jets.size();
1990 0 : vector<EtaPhi> points(n); // recall EtaPhi is just a typedef'd pair<double>
1991 0 : for (int i = 0; i < n; i++) {
1992 0 : points[i] = EtaPhi(_jets[i].rap(),_jets[i].phi_02pi());
1993 0 : points[i].sanitize(); // make sure things are in the right range
1994 : }
1995 0 : auto_ptr<DynamicNearestNeighbours> DNN;
1996 : #ifndef __FJCORE_DROP_CGAL // strategy = NlnN* are not supported if we drop CGAL...
1997 : bool verbose = false;
1998 : bool ignore_nearest_is_mirror = (_Rparam < twopi);
1999 : if (_strategy == NlnN4pi) {
2000 : DNN.reset(new Dnn4piCylinder(points,verbose));
2001 : } else if (_strategy == NlnN3pi) {
2002 : DNN.reset(new Dnn3piCylinder(points,ignore_nearest_is_mirror,verbose));
2003 : } else if (_strategy == NlnN) {
2004 : DNN.reset(new Dnn2piCylinder(points,ignore_nearest_is_mirror,verbose));
2005 : } else
2006 : #else
2007 0 : if (_strategy == NlnN4pi || _strategy == NlnN3pi || _strategy == NlnN) {
2008 0 : ostringstream err;
2009 0 : err << "ERROR: Requested strategy "<<strategy_string()<<" but it is not"<<endl;
2010 0 : err << " supported because FastJet was compiled without CGAL"<<endl;
2011 0 : throw Error(err.str());
2012 0 : }
2013 : #endif // __FJCORE_DROP_CGAL
2014 : {
2015 0 : ostringstream err;
2016 0 : err << "ERROR: Unrecognized value for strategy: "<<_strategy<<endl;
2017 0 : assert(false);
2018 : throw Error(err.str());
2019 0 : }
2020 : DistMap DijMap;
2021 : for (int ii = 0; ii < n; ii++) {
2022 : _add_ktdistance_to_map(ii, DijMap, DNN.get());
2023 : }
2024 : for (int i=0;i<n;i++) {
2025 : TwoVertices SmallestDijPair;
2026 : int jet_i, jet_j;
2027 : double SmallestDij;
2028 : bool Valid2;
2029 : bool recombine_with_beam;
2030 : do {
2031 : SmallestDij = DijMap.begin()->first;
2032 : SmallestDijPair = DijMap.begin()->second;
2033 : jet_i = SmallestDijPair.first;
2034 : jet_j = SmallestDijPair.second;
2035 : DijMap.erase(DijMap.begin());
2036 : recombine_with_beam = (jet_j == BeamJet);
2037 : if (!recombine_with_beam) {Valid2 = DNN->Valid(jet_j);}
2038 : else {Valid2 = true;}
2039 : } while ( !DNN->Valid(jet_i) || !Valid2);
2040 : if (! recombine_with_beam) {
2041 : int nn; // will be index of new jet
2042 : _do_ij_recombination_step(jet_i, jet_j, SmallestDij, nn);
2043 : EtaPhi newpoint(_jets[nn].rap(), _jets[nn].phi_02pi());
2044 : newpoint.sanitize(); // make sure it is in correct range
2045 : points.push_back(newpoint);
2046 : } else {
2047 : _do_iB_recombination_step(jet_i, SmallestDij);
2048 : }
2049 : if (i == n-1) {break;}
2050 : vector<int> updated_neighbours;
2051 : if (! recombine_with_beam) {
2052 : int point3;
2053 : DNN->RemoveCombinedAddCombination(jet_i, jet_j,
2054 : points[points.size()-1], point3,
2055 : updated_neighbours);
2056 : if (static_cast<unsigned int> (point3) != points.size()-1) {
2057 : throw Error("INTERNAL ERROR: point3 != points.size()-1");}
2058 : } else {
2059 : DNN->RemovePoint(jet_i, updated_neighbours);
2060 : }
2061 : vector<int>::iterator it = updated_neighbours.begin();
2062 : for (; it != updated_neighbours.end(); ++it) {
2063 : int ii = *it;
2064 : _add_ktdistance_to_map(ii, DijMap, DNN.get());
2065 : }
2066 : } // end clustering loop
2067 0 : }
2068 : void ClusterSequence::_add_ktdistance_to_map(
2069 : const int & ii,
2070 : DistMap & DijMap,
2071 : const DynamicNearestNeighbours * DNN) {
2072 0 : double yiB = jet_scale_for_algorithm(_jets[ii]);
2073 0 : if (yiB == 0.0) {
2074 0 : DijMap.insert(DijEntry(yiB, TwoVertices(ii,-1)));
2075 0 : } else {
2076 0 : double DeltaR2 = DNN->NearestNeighbourDistance(ii) * _invR2;
2077 0 : if (DeltaR2 > 1.0) {
2078 0 : DijMap.insert(DijEntry(yiB, TwoVertices(ii,-1)));
2079 0 : } else {
2080 0 : double kt2i = jet_scale_for_algorithm(_jets[ii]);
2081 0 : int jj = DNN->NearestNeighbourIndex(ii);
2082 0 : if (kt2i <= jet_scale_for_algorithm(_jets[jj])) {
2083 0 : double dij = DeltaR2 * kt2i;
2084 0 : DijMap.insert(DijEntry(dij, TwoVertices(ii,jj)));
2085 0 : }
2086 0 : }
2087 : }
2088 0 : }
2089 : FJCORE_END_NAMESPACE
2090 : #include<iostream>
2091 : #include<cmath>
2092 : #include <cstdlib>
2093 : #include<cassert>
2094 : FJCORE_BEGIN_NAMESPACE // defined in fastjet/internal/base.hh
2095 : using namespace std;
2096 : void ClusterSequence::_really_dumb_cluster () {
2097 0 : vector<PseudoJet *> jetsp(_jets.size());
2098 0 : vector<int> indices(_jets.size());
2099 0 : for (size_t i = 0; i<_jets.size(); i++) {
2100 0 : jetsp[i] = & _jets[i];
2101 0 : indices[i] = i;
2102 : }
2103 0 : for (int n = jetsp.size(); n > 0; n--) {
2104 : int ii, jj;
2105 0 : double ymin = jet_scale_for_algorithm(*(jetsp[0]));
2106 : ii = 0; jj = -2;
2107 0 : for (int i = 0; i < n; i++) {
2108 0 : double yiB = jet_scale_for_algorithm(*(jetsp[i]));
2109 0 : if (yiB < ymin) {
2110 0 : ymin = yiB; ii = i; jj = -2;}
2111 : }
2112 0 : for (int i = 0; i < n-1; i++) {
2113 0 : for (int j = i+1; j < n; j++) {
2114 : //double y = jetsp[i]->kt_distance(*jetsp[j])*_invR2;
2115 0 : double y = min(jet_scale_for_algorithm(*(jetsp[i])),
2116 0 : jet_scale_for_algorithm(*(jetsp[j])))
2117 0 : * jetsp[i]->plain_distance(*jetsp[j])*_invR2;
2118 0 : if (y < ymin) {ymin = y; ii = i; jj = j;}
2119 : }
2120 : }
2121 0 : int newn = 2*jetsp.size() - n;
2122 0 : if (jj >= 0) {
2123 0 : int nn; // new jet index
2124 0 : _do_ij_recombination_step(jetsp[ii]-&_jets[0],
2125 0 : jetsp[jj]-&_jets[0], ymin, nn);
2126 0 : jetsp[ii] = &_jets[nn];
2127 0 : jetsp[jj] = jetsp[n-1];
2128 0 : indices[ii] = newn;
2129 0 : indices[jj] = indices[n-1];
2130 0 : } else {
2131 0 : _do_iB_recombination_step(jetsp[ii]-&_jets[0], ymin);
2132 0 : jetsp[ii] = jetsp[n-1];
2133 0 : indices[ii] = indices[n-1];
2134 : }
2135 0 : }
2136 0 : }
2137 : FJCORE_END_NAMESPACE
2138 : #include<iostream>
2139 : FJCORE_BEGIN_NAMESPACE // defined in fastjet/internal/base.hh
2140 : using namespace std;
2141 : template<> inline void ClusterSequence::_bj_set_jetinfo(
2142 : EEBriefJet * const jetA, const int _jets_index) const {
2143 0 : double E = _jets[_jets_index].E();
2144 0 : double scale = E*E; // the default energy scale for the kt alg
2145 0 : double p = jet_def().extra_param(); // in case we're ee_genkt
2146 0 : switch (_jet_algorithm) {
2147 : case ee_kt_algorithm:
2148 0 : assert(_Rparam > 2.0); // force this to be true! [not best place, but works]
2149 : break;
2150 : case ee_genkt_algorithm:
2151 0 : if (p <= 0 && scale < 1e-300) scale = 1e-300; // same dodgy safety as genkt
2152 0 : scale = pow(scale,p);
2153 0 : break;
2154 : default:
2155 0 : throw Error("Unrecognised jet algorithm");
2156 : }
2157 0 : jetA->kt2 = scale; // "kt2" might one day be renamed as "scale" or some such
2158 0 : double norm = _jets[_jets_index].modp2();
2159 0 : if (norm > 0) {
2160 0 : norm = 1.0/sqrt(norm);
2161 0 : jetA->nx = norm * _jets[_jets_index].px();
2162 0 : jetA->ny = norm * _jets[_jets_index].py();
2163 0 : jetA->nz = norm * _jets[_jets_index].pz();
2164 0 : } else {
2165 0 : jetA->nx = 0.0;
2166 0 : jetA->ny = 0.0;
2167 0 : jetA->nz = 1.0;
2168 : }
2169 0 : jetA->_jets_index = _jets_index;
2170 0 : jetA->NN_dist = _R2;
2171 0 : jetA->NN = NULL;
2172 0 : }
2173 : template<> double ClusterSequence::_bj_dist(
2174 : const EEBriefJet * const jeta,
2175 : const EEBriefJet * const jetb) const {
2176 : double dist = 1.0
2177 0 : - jeta->nx*jetb->nx
2178 0 : - jeta->ny*jetb->ny
2179 0 : - jeta->nz*jetb->nz;
2180 0 : dist *= 2; // distance is _2_*min(Ei^2,Ej^2)*(1-cos theta)
2181 0 : return dist;
2182 : }
2183 : void ClusterSequence::_simple_N2_cluster_BriefJet() {
2184 0 : _simple_N2_cluster<BriefJet>();
2185 0 : }
2186 : void ClusterSequence::_simple_N2_cluster_EEBriefJet() {
2187 0 : _simple_N2_cluster<EEBriefJet>();
2188 0 : }
2189 : FJCORE_END_NAMESPACE
2190 : #include <iostream>
2191 : FJCORE_BEGIN_NAMESPACE // defined in fastjet/internal/base.hh
2192 : using namespace std;
2193 0 : ClusterSequenceStructure::~ClusterSequenceStructure(){
2194 0 : if (_associated_cs != NULL
2195 0 : && _associated_cs->will_delete_self_when_unused()) {
2196 0 : _associated_cs->signal_imminent_self_deletion();
2197 0 : delete _associated_cs;
2198 : }
2199 0 : }
2200 : bool ClusterSequenceStructure::has_valid_cluster_sequence() const{
2201 0 : return (_associated_cs != NULL);
2202 : }
2203 : const ClusterSequence* ClusterSequenceStructure::associated_cluster_sequence() const{
2204 0 : return _associated_cs;
2205 : }
2206 : const ClusterSequence * ClusterSequenceStructure::validated_cs() const {
2207 0 : if (!_associated_cs)
2208 0 : throw Error("you requested information about the internal structure of a jet, but its associated ClusterSequence has gone out of scope.");
2209 0 : return _associated_cs;
2210 0 : }
2211 : bool ClusterSequenceStructure::has_partner(const PseudoJet &reference, PseudoJet &partner) const{
2212 0 : return validated_cs()->has_partner(reference, partner);
2213 : }
2214 : bool ClusterSequenceStructure::has_child(const PseudoJet &reference, PseudoJet &child) const{
2215 0 : return validated_cs()->has_child(reference, child);
2216 : }
2217 : bool ClusterSequenceStructure::has_parents(const PseudoJet &reference, PseudoJet &parent1, PseudoJet &parent2) const{
2218 0 : return validated_cs()->has_parents(reference, parent1, parent2);
2219 : }
2220 : bool ClusterSequenceStructure::object_in_jet(const PseudoJet &reference, const PseudoJet &jet) const{
2221 0 : if ((!has_associated_cluster_sequence()) || (!jet.has_associated_cluster_sequence()))
2222 0 : 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.");
2223 0 : if (reference.associated_cluster_sequence() != jet.associated_cluster_sequence()) return false;
2224 0 : return validated_cs()->object_in_jet(reference, jet);
2225 0 : }
2226 : bool ClusterSequenceStructure::has_constituents() const{
2227 0 : if (!has_associated_cluster_sequence())
2228 0 : 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.");
2229 0 : return true;
2230 0 : }
2231 : vector<PseudoJet> ClusterSequenceStructure::constituents(const PseudoJet &reference) const{
2232 0 : return validated_cs()->constituents(reference);
2233 : }
2234 : bool ClusterSequenceStructure::has_exclusive_subjets() const{
2235 0 : if (!has_associated_cluster_sequence())
2236 0 : 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.");
2237 0 : return true;
2238 0 : }
2239 : std::vector<PseudoJet> ClusterSequenceStructure::exclusive_subjets (const PseudoJet &reference, const double & dcut) const {
2240 0 : return validated_cs()->exclusive_subjets(reference, dcut);
2241 : }
2242 : int ClusterSequenceStructure::n_exclusive_subjets(const PseudoJet &reference, const double & dcut) const {
2243 0 : return validated_cs()->n_exclusive_subjets(reference, dcut);
2244 : }
2245 : std::vector<PseudoJet> ClusterSequenceStructure::exclusive_subjets_up_to (const PseudoJet &reference, int nsub) const {
2246 0 : return validated_cs()->exclusive_subjets_up_to(reference, nsub);
2247 : }
2248 : double ClusterSequenceStructure::exclusive_subdmerge(const PseudoJet &reference, int nsub) const {
2249 0 : return validated_cs()->exclusive_subdmerge(reference, nsub);
2250 : }
2251 : double ClusterSequenceStructure::exclusive_subdmerge_max(const PseudoJet &reference, int nsub) const {
2252 0 : return validated_cs()->exclusive_subdmerge_max(reference, nsub);
2253 : }
2254 : bool ClusterSequenceStructure::has_pieces(const PseudoJet &reference) const{
2255 0 : PseudoJet dummy1, dummy2;
2256 0 : return has_parents(reference, dummy1, dummy2);
2257 0 : }
2258 : vector<PseudoJet> ClusterSequenceStructure::pieces(const PseudoJet &reference) const{
2259 0 : PseudoJet j1, j2;
2260 0 : vector<PseudoJet> res;
2261 0 : if (has_parents(reference, j1, j2)){
2262 0 : res.push_back(j1);
2263 0 : res.push_back(j2);
2264 : }
2265 : return res;
2266 0 : }
2267 : FJCORE_END_NAMESPACE
2268 : #include<iostream>
2269 : #include<vector>
2270 : #include<cmath>
2271 : #include<algorithm>
2272 : FJCORE_BEGIN_NAMESPACE // defined in fastjet/internal/base.hh
2273 : using namespace std;
2274 : void ClusterSequence::_bj_remove_from_tiles(TiledJet * const jet) {
2275 0 : Tile * tile = & _tiles[jet->tile_index];
2276 0 : if (jet->previous == NULL) {
2277 0 : tile->head = jet->next;
2278 0 : } else {
2279 0 : jet->previous->next = jet->next;
2280 : }
2281 0 : if (jet->next != NULL) {
2282 0 : jet->next->previous = jet->previous;
2283 0 : }
2284 0 : }
2285 : void ClusterSequence::_initialise_tiles() {
2286 0 : double default_size = max(0.1,_Rparam);
2287 0 : _tile_size_eta = default_size;
2288 0 : _n_tiles_phi = max(3,int(floor(twopi/default_size)));
2289 0 : _tile_size_phi = twopi / _n_tiles_phi; // >= _Rparam and fits in 2pi
2290 0 : _tiles_eta_min = 0.0;
2291 0 : _tiles_eta_max = 0.0;
2292 : const double maxrap = 7.0;
2293 0 : for(unsigned int i = 0; i < _jets.size(); i++) {
2294 0 : double eta = _jets[i].rap();
2295 0 : if (abs(eta) < maxrap) {
2296 0 : if (eta < _tiles_eta_min) {_tiles_eta_min = eta;}
2297 0 : if (eta > _tiles_eta_max) {_tiles_eta_max = eta;}
2298 : }
2299 : }
2300 0 : _tiles_ieta_min = int(floor(_tiles_eta_min/_tile_size_eta));
2301 0 : _tiles_ieta_max = int(floor( _tiles_eta_max/_tile_size_eta));
2302 0 : _tiles_eta_min = _tiles_ieta_min * _tile_size_eta;
2303 0 : _tiles_eta_max = _tiles_ieta_max * _tile_size_eta;
2304 0 : _tiles.resize((_tiles_ieta_max-_tiles_ieta_min+1)*_n_tiles_phi);
2305 0 : for (int ieta = _tiles_ieta_min; ieta <= _tiles_ieta_max; ieta++) {
2306 0 : for (int iphi = 0; iphi < _n_tiles_phi; iphi++) {
2307 0 : Tile * tile = & _tiles[_tile_index(ieta,iphi)];
2308 0 : tile->head = NULL; // first element of tiles points to itself
2309 0 : tile->begin_tiles[0] = tile;
2310 : Tile ** pptile = & (tile->begin_tiles[0]);
2311 0 : pptile++;
2312 0 : tile->surrounding_tiles = pptile;
2313 0 : if (ieta > _tiles_ieta_min) {
2314 : // with the itile subroutine, we can safely run tiles from
2315 : // idphi=-1 to idphi=+1, because it takes care of
2316 : // negative and positive boundaries
2317 0 : for (int idphi = -1; idphi <=+1; idphi++) {
2318 0 : *pptile = & _tiles[_tile_index(ieta-1,iphi+idphi)];
2319 0 : pptile++;
2320 : }
2321 0 : }
2322 0 : *pptile = & _tiles[_tile_index(ieta,iphi-1)];
2323 0 : pptile++;
2324 0 : tile->RH_tiles = pptile;
2325 0 : *pptile = & _tiles[_tile_index(ieta,iphi+1)];
2326 0 : pptile++;
2327 0 : if (ieta < _tiles_ieta_max) {
2328 0 : for (int idphi = -1; idphi <= +1; idphi++) {
2329 0 : *pptile = & _tiles[_tile_index(ieta+1,iphi+idphi)];
2330 0 : pptile++;
2331 : }
2332 0 : }
2333 0 : tile->end_tiles = pptile;
2334 0 : tile->tagged = false;
2335 : }
2336 : }
2337 0 : }
2338 : int ClusterSequence::_tile_index(const double & eta, const double & phi) const {
2339 : int ieta, iphi;
2340 0 : if (eta <= _tiles_eta_min) {ieta = 0;}
2341 0 : else if (eta >= _tiles_eta_max) {ieta = _tiles_ieta_max-_tiles_ieta_min;}
2342 : else {
2343 0 : ieta = int(((eta - _tiles_eta_min) / _tile_size_eta));
2344 0 : if (ieta > _tiles_ieta_max-_tiles_ieta_min) {
2345 0 : ieta = _tiles_ieta_max-_tiles_ieta_min;}
2346 : }
2347 0 : iphi = int((phi+twopi)/_tile_size_phi) % _n_tiles_phi;
2348 0 : return (iphi + ieta * _n_tiles_phi);
2349 : }
2350 : inline void ClusterSequence::_tj_set_jetinfo( TiledJet * const jet,
2351 : const int _jets_index) {
2352 0 : _bj_set_jetinfo<>(jet, _jets_index);
2353 0 : jet->tile_index = _tile_index(jet->eta, jet->phi);
2354 0 : Tile * tile = &_tiles[jet->tile_index];
2355 0 : jet->previous = NULL;
2356 0 : jet->next = tile->head;
2357 0 : if (jet->next != NULL) {jet->next->previous = jet;}
2358 0 : tile->head = jet;
2359 0 : }
2360 : void ClusterSequence::_print_tiles(TiledJet * briefjets ) const {
2361 0 : for (vector<Tile>::const_iterator tile = _tiles.begin();
2362 0 : tile < _tiles.end(); tile++) {
2363 0 : cout << "Tile " << tile - _tiles.begin()<<" = ";
2364 0 : vector<int> list;
2365 0 : for (TiledJet * jetI = tile->head; jetI != NULL; jetI = jetI->next) {
2366 0 : list.push_back(jetI-briefjets);
2367 : }
2368 0 : sort(list.begin(),list.end());
2369 0 : for (unsigned int i = 0; i < list.size(); i++) {cout <<" "<<list[i];}
2370 0 : cout <<"\n";
2371 0 : }
2372 0 : }
2373 : void ClusterSequence::_add_neighbours_to_tile_union(const int tile_index,
2374 : vector<int> & tile_union, int & n_near_tiles) const {
2375 0 : for (Tile * const * near_tile = _tiles[tile_index].begin_tiles;
2376 0 : near_tile != _tiles[tile_index].end_tiles; near_tile++){
2377 0 : tile_union[n_near_tiles] = *near_tile - & _tiles[0];
2378 0 : n_near_tiles++;
2379 : }
2380 0 : }
2381 : inline void ClusterSequence::_add_untagged_neighbours_to_tile_union(
2382 : const int tile_index,
2383 : vector<int> & tile_union, int & n_near_tiles) {
2384 0 : for (Tile ** near_tile = _tiles[tile_index].begin_tiles;
2385 0 : near_tile != _tiles[tile_index].end_tiles; near_tile++){
2386 0 : if (! (*near_tile)->tagged) {
2387 0 : (*near_tile)->tagged = true;
2388 0 : tile_union[n_near_tiles] = *near_tile - & _tiles[0];
2389 0 : n_near_tiles++;
2390 0 : }
2391 : }
2392 0 : }
2393 : void ClusterSequence::_tiled_N2_cluster() {
2394 0 : _initialise_tiles();
2395 0 : int n = _jets.size();
2396 0 : TiledJet * briefjets = new TiledJet[n];
2397 0 : TiledJet * jetA = briefjets, * jetB;
2398 0 : TiledJet oldB;
2399 : oldB.tile_index=0; // prevents a gcc warning
2400 0 : vector<int> tile_union(3*n_tile_neighbours);
2401 0 : for (int i = 0; i< n; i++) {
2402 0 : _tj_set_jetinfo(jetA, i);
2403 0 : jetA++; // move on to next entry of briefjets
2404 : }
2405 0 : TiledJet * tail = jetA; // a semaphore for the end of briefjets
2406 : TiledJet * head = briefjets; // a nicer way of naming start
2407 0 : vector<Tile>::const_iterator tile;
2408 0 : for (tile = _tiles.begin(); tile != _tiles.end(); tile++) {
2409 0 : for (jetA = tile->head; jetA != NULL; jetA = jetA->next) {
2410 0 : for (jetB = tile->head; jetB != jetA; jetB = jetB->next) {
2411 0 : double dist = _bj_dist(jetA,jetB);
2412 0 : if (dist < jetA->NN_dist) {jetA->NN_dist = dist; jetA->NN = jetB;}
2413 0 : if (dist < jetB->NN_dist) {jetB->NN_dist = dist; jetB->NN = jetA;}
2414 : }
2415 : }
2416 0 : for (Tile ** RTile = tile->RH_tiles; RTile != tile->end_tiles; RTile++) {
2417 0 : for (jetA = tile->head; jetA != NULL; jetA = jetA->next) {
2418 0 : for (jetB = (*RTile)->head; jetB != NULL; jetB = jetB->next) {
2419 0 : double dist = _bj_dist(jetA,jetB);
2420 0 : if (dist < jetA->NN_dist) {jetA->NN_dist = dist; jetA->NN = jetB;}
2421 0 : if (dist < jetB->NN_dist) {jetB->NN_dist = dist; jetB->NN = jetA;}
2422 : }
2423 : }
2424 : }
2425 : }
2426 0 : double * diJ = new double[n];
2427 0 : jetA = head;
2428 0 : for (int i = 0; i < n; i++) {
2429 0 : diJ[i] = _bj_diJ(jetA);
2430 0 : jetA++; // have jetA follow i
2431 : }
2432 0 : int history_location = n-1;
2433 0 : while (tail != head) {
2434 0 : double diJ_min = diJ[0];
2435 : int diJ_min_jet = 0;
2436 0 : for (int i = 1; i < n; i++) {
2437 0 : if (diJ[i] < diJ_min) {diJ_min_jet = i; diJ_min = diJ[i];}
2438 : }
2439 0 : history_location++;
2440 0 : jetA = & briefjets[diJ_min_jet];
2441 0 : jetB = jetA->NN;
2442 0 : diJ_min *= _invR2;
2443 0 : if (jetB != NULL) {
2444 0 : if (jetA < jetB) {std::swap(jetA,jetB);}
2445 0 : int nn; // new jet index
2446 0 : _do_ij_recombination_step(jetA->_jets_index, jetB->_jets_index, diJ_min, nn);
2447 0 : _bj_remove_from_tiles(jetA);
2448 0 : oldB = * jetB; // take a copy because we will need it...
2449 0 : _bj_remove_from_tiles(jetB);
2450 0 : _tj_set_jetinfo(jetB, nn); // also registers the jet in the tiling
2451 0 : } else {
2452 0 : _do_iB_recombination_step(jetA->_jets_index, diJ_min);
2453 0 : _bj_remove_from_tiles(jetA);
2454 : }
2455 0 : int n_near_tiles = 0;
2456 0 : _add_neighbours_to_tile_union(jetA->tile_index, tile_union, n_near_tiles);
2457 0 : if (jetB != NULL) {
2458 : bool sort_it = false;
2459 0 : if (jetB->tile_index != jetA->tile_index) {
2460 : sort_it = true;
2461 0 : _add_neighbours_to_tile_union(jetB->tile_index,tile_union,n_near_tiles);
2462 : }
2463 0 : if (oldB.tile_index != jetA->tile_index &&
2464 0 : oldB.tile_index != jetB->tile_index) {
2465 : sort_it = true;
2466 0 : _add_neighbours_to_tile_union(oldB.tile_index,tile_union,n_near_tiles);
2467 : }
2468 0 : if (sort_it) {
2469 : // sort the tiles before then compressing the list
2470 0 : sort(tile_union.begin(), tile_union.begin()+n_near_tiles);
2471 : // and now condense the list
2472 : int nnn = 1;
2473 0 : for (int i = 1; i < n_near_tiles; i++) {
2474 0 : if (tile_union[i] != tile_union[nnn-1]) {
2475 0 : tile_union[nnn] = tile_union[i];
2476 0 : nnn++;
2477 0 : }
2478 : }
2479 0 : n_near_tiles = nnn;
2480 0 : }
2481 0 : }
2482 0 : tail--; n--;
2483 0 : if (jetA == tail) {
2484 : } else {
2485 0 : *jetA = *tail;
2486 0 : diJ[jetA - head] = diJ[tail-head];
2487 0 : if (jetA->previous == NULL) {
2488 0 : _tiles[jetA->tile_index].head = jetA;
2489 0 : } else {
2490 0 : jetA->previous->next = jetA;
2491 : }
2492 0 : if (jetA->next != NULL) {jetA->next->previous = jetA;}
2493 : }
2494 0 : for (int itile = 0; itile < n_near_tiles; itile++) {
2495 0 : Tile * tile_ptr = &_tiles[tile_union[itile]];
2496 0 : for (TiledJet * jetI = tile_ptr->head; jetI != NULL; jetI = jetI->next) {
2497 : // see if jetI had jetA or jetB as a NN -- if so recalculate the NN
2498 0 : if (jetI->NN == jetA || (jetI->NN == jetB && jetB != NULL)) {
2499 0 : jetI->NN_dist = _R2;
2500 0 : jetI->NN = NULL;
2501 : // now go over tiles that are neighbours of I (include own tile)
2502 0 : for (Tile ** near_tile = tile_ptr->begin_tiles;
2503 0 : near_tile != tile_ptr->end_tiles; near_tile++) {
2504 : // and then over the contents of that tile
2505 0 : for (TiledJet * jetJ = (*near_tile)->head;
2506 0 : jetJ != NULL; jetJ = jetJ->next) {
2507 0 : double dist = _bj_dist(jetI,jetJ);
2508 0 : if (dist < jetI->NN_dist && jetJ != jetI) {
2509 0 : jetI->NN_dist = dist; jetI->NN = jetJ;
2510 0 : }
2511 : }
2512 : }
2513 0 : diJ[jetI-head] = _bj_diJ(jetI); // update diJ
2514 0 : }
2515 : // check whether new jetB is closer than jetI's current NN and
2516 : // if need to update things
2517 0 : if (jetB != NULL) {
2518 0 : double dist = _bj_dist(jetI,jetB);
2519 0 : if (dist < jetI->NN_dist) {
2520 0 : if (jetI != jetB) {
2521 0 : jetI->NN_dist = dist;
2522 0 : jetI->NN = jetB;
2523 0 : diJ[jetI-head] = _bj_diJ(jetI); // update diJ...
2524 0 : }
2525 : }
2526 0 : if (dist < jetB->NN_dist) {
2527 0 : if (jetI != jetB) {
2528 0 : jetB->NN_dist = dist;
2529 0 : jetB->NN = jetI;}
2530 : }
2531 0 : }
2532 : }
2533 : }
2534 0 : if (jetB != NULL) {diJ[jetB-head] = _bj_diJ(jetB);}
2535 0 : for (Tile ** near_tile = _tiles[tail->tile_index].begin_tiles;
2536 0 : near_tile!= _tiles[tail->tile_index].end_tiles; near_tile++){
2537 0 : for (TiledJet * jetJ = (*near_tile)->head;
2538 0 : jetJ != NULL; jetJ = jetJ->next) {
2539 0 : if (jetJ->NN == tail) {jetJ->NN = jetA;}
2540 : }
2541 : }
2542 0 : if (jetB != NULL) {diJ[jetB-head] = _bj_diJ(jetB);}
2543 0 : }
2544 0 : delete[] diJ;
2545 0 : delete[] briefjets;
2546 0 : }
2547 : void ClusterSequence::_faster_tiled_N2_cluster() {
2548 0 : _initialise_tiles();
2549 0 : int n = _jets.size();
2550 0 : TiledJet * briefjets = new TiledJet[n];
2551 0 : TiledJet * jetA = briefjets, * jetB;
2552 0 : TiledJet oldB;
2553 : oldB.tile_index=0; // prevents a gcc warning
2554 0 : vector<int> tile_union(3*n_tile_neighbours);
2555 0 : for (int i = 0; i< n; i++) {
2556 0 : _tj_set_jetinfo(jetA, i);
2557 0 : jetA++; // move on to next entry of briefjets
2558 : }
2559 : TiledJet * head = briefjets; // a nicer way of naming start
2560 0 : vector<Tile>::const_iterator tile;
2561 0 : for (tile = _tiles.begin(); tile != _tiles.end(); tile++) {
2562 0 : for (jetA = tile->head; jetA != NULL; jetA = jetA->next) {
2563 0 : for (jetB = tile->head; jetB != jetA; jetB = jetB->next) {
2564 0 : double dist = _bj_dist(jetA,jetB);
2565 0 : if (dist < jetA->NN_dist) {jetA->NN_dist = dist; jetA->NN = jetB;}
2566 0 : if (dist < jetB->NN_dist) {jetB->NN_dist = dist; jetB->NN = jetA;}
2567 : }
2568 : }
2569 0 : for (Tile ** RTile = tile->RH_tiles; RTile != tile->end_tiles; RTile++) {
2570 0 : for (jetA = tile->head; jetA != NULL; jetA = jetA->next) {
2571 0 : for (jetB = (*RTile)->head; jetB != NULL; jetB = jetB->next) {
2572 0 : double dist = _bj_dist(jetA,jetB);
2573 0 : if (dist < jetA->NN_dist) {jetA->NN_dist = dist; jetA->NN = jetB;}
2574 0 : if (dist < jetB->NN_dist) {jetB->NN_dist = dist; jetB->NN = jetA;}
2575 : }
2576 : }
2577 : }
2578 : }
2579 : struct diJ_plus_link {
2580 : double diJ; // the distance
2581 : TiledJet * jet; // the jet (i) for which we've found this distance
2582 : };
2583 0 : diJ_plus_link * diJ = new diJ_plus_link[n];
2584 0 : jetA = head;
2585 0 : for (int i = 0; i < n; i++) {
2586 0 : diJ[i].diJ = _bj_diJ(jetA); // kt distance * R^2
2587 0 : diJ[i].jet = jetA; // our compact diJ table will not be in
2588 0 : jetA->diJ_posn = i; // one-to-one corresp. with non-compact jets,
2589 0 : jetA++; // have jetA follow i
2590 : }
2591 0 : int history_location = n-1;
2592 0 : while (n > 0) {
2593 : diJ_plus_link * best, *stop; // pointers a bit faster than indices
2594 0 : double diJ_min = diJ[0].diJ; // initialise the best one here.
2595 : best = diJ; // and here
2596 0 : stop = diJ+n;
2597 0 : for (diJ_plus_link * here = diJ+1; here != stop; here++) {
2598 0 : if (here->diJ < diJ_min) {best = here; diJ_min = here->diJ;}
2599 : }
2600 0 : history_location++;
2601 0 : jetA = best->jet;
2602 0 : jetB = jetA->NN;
2603 0 : diJ_min *= _invR2;
2604 0 : if (jetB != NULL) {
2605 0 : if (jetA < jetB) {std::swap(jetA,jetB);}
2606 0 : int nn; // new jet index
2607 0 : _do_ij_recombination_step(jetA->_jets_index, jetB->_jets_index, diJ_min, nn);
2608 0 : _bj_remove_from_tiles(jetA);
2609 0 : oldB = * jetB; // take a copy because we will need it...
2610 0 : _bj_remove_from_tiles(jetB);
2611 0 : _tj_set_jetinfo(jetB, nn); // cause jetB to become _jets[nn]
2612 0 : } else {
2613 0 : _do_iB_recombination_step(jetA->_jets_index, diJ_min);
2614 0 : _bj_remove_from_tiles(jetA);
2615 : }
2616 0 : int n_near_tiles = 0;
2617 0 : _add_untagged_neighbours_to_tile_union(jetA->tile_index,
2618 : tile_union, n_near_tiles);
2619 0 : if (jetB != NULL) {
2620 0 : if (jetB->tile_index != jetA->tile_index) {
2621 0 : _add_untagged_neighbours_to_tile_union(jetB->tile_index,
2622 : tile_union,n_near_tiles);
2623 : }
2624 0 : if (oldB.tile_index != jetA->tile_index &&
2625 0 : oldB.tile_index != jetB->tile_index) {
2626 0 : _add_untagged_neighbours_to_tile_union(oldB.tile_index,
2627 : tile_union,n_near_tiles);
2628 : }
2629 : }
2630 0 : n--;
2631 0 : diJ[n].jet->diJ_posn = jetA->diJ_posn;
2632 0 : diJ[jetA->diJ_posn] = diJ[n];
2633 0 : for (int itile = 0; itile < n_near_tiles; itile++) {
2634 0 : Tile * tile_ptr = &_tiles[tile_union[itile]];
2635 0 : tile_ptr->tagged = false; // reset tag, since we're done with unions
2636 0 : for (TiledJet * jetI = tile_ptr->head; jetI != NULL; jetI = jetI->next) {
2637 : // see if jetI had jetA or jetB as a NN -- if so recalculate the NN
2638 0 : if (jetI->NN == jetA || (jetI->NN == jetB && jetB != NULL)) {
2639 0 : jetI->NN_dist = _R2;
2640 0 : jetI->NN = NULL;
2641 : // now go over tiles that are neighbours of I (include own tile)
2642 0 : for (Tile ** near_tile = tile_ptr->begin_tiles;
2643 0 : near_tile != tile_ptr->end_tiles; near_tile++) {
2644 : // and then over the contents of that tile
2645 0 : for (TiledJet * jetJ = (*near_tile)->head;
2646 0 : jetJ != NULL; jetJ = jetJ->next) {
2647 0 : double dist = _bj_dist(jetI,jetJ);
2648 0 : if (dist < jetI->NN_dist && jetJ != jetI) {
2649 0 : jetI->NN_dist = dist; jetI->NN = jetJ;
2650 0 : }
2651 : }
2652 : }
2653 0 : diJ[jetI->diJ_posn].diJ = _bj_diJ(jetI); // update diJ kt-dist
2654 0 : }
2655 : // check whether new jetB is closer than jetI's current NN and
2656 : // if jetI is closer than jetB's current (evolving) nearest
2657 : // neighbour. Where relevant update things
2658 0 : if (jetB != NULL) {
2659 0 : double dist = _bj_dist(jetI,jetB);
2660 0 : if (dist < jetI->NN_dist) {
2661 0 : if (jetI != jetB) {
2662 0 : jetI->NN_dist = dist;
2663 0 : jetI->NN = jetB;
2664 0 : diJ[jetI->diJ_posn].diJ = _bj_diJ(jetI); // update diJ...
2665 0 : }
2666 : }
2667 0 : if (dist < jetB->NN_dist) {
2668 0 : if (jetI != jetB) {
2669 0 : jetB->NN_dist = dist;
2670 0 : jetB->NN = jetI;}
2671 : }
2672 0 : }
2673 : }
2674 : }
2675 0 : if (jetB != NULL) {diJ[jetB->diJ_posn].diJ = _bj_diJ(jetB);}
2676 0 : }
2677 0 : delete[] diJ;
2678 0 : delete[] briefjets;
2679 0 : }
2680 : void ClusterSequence::_minheap_faster_tiled_N2_cluster() {
2681 0 : _initialise_tiles();
2682 0 : int n = _jets.size();
2683 0 : TiledJet * briefjets = new TiledJet[n];
2684 0 : TiledJet * jetA = briefjets, * jetB;
2685 0 : TiledJet oldB;
2686 : oldB.tile_index=0; // prevents a gcc warning
2687 0 : vector<int> tile_union(3*n_tile_neighbours);
2688 0 : for (int i = 0; i< n; i++) {
2689 0 : _tj_set_jetinfo(jetA, i);
2690 0 : jetA++; // move on to next entry of briefjets
2691 : }
2692 : TiledJet * head = briefjets; // a nicer way of naming start
2693 0 : vector<Tile>::const_iterator tile;
2694 0 : for (tile = _tiles.begin(); tile != _tiles.end(); tile++) {
2695 0 : for (jetA = tile->head; jetA != NULL; jetA = jetA->next) {
2696 0 : for (jetB = tile->head; jetB != jetA; jetB = jetB->next) {
2697 0 : double dist = _bj_dist(jetA,jetB);
2698 0 : if (dist < jetA->NN_dist) {jetA->NN_dist = dist; jetA->NN = jetB;}
2699 0 : if (dist < jetB->NN_dist) {jetB->NN_dist = dist; jetB->NN = jetA;}
2700 : }
2701 : }
2702 0 : for (Tile ** RTile = tile->RH_tiles; RTile != tile->end_tiles; RTile++) {
2703 0 : for (jetA = tile->head; jetA != NULL; jetA = jetA->next) {
2704 0 : for (jetB = (*RTile)->head; jetB != NULL; jetB = jetB->next) {
2705 0 : double dist = _bj_dist(jetA,jetB);
2706 0 : if (dist < jetA->NN_dist) {jetA->NN_dist = dist; jetA->NN = jetB;}
2707 0 : if (dist < jetB->NN_dist) {jetB->NN_dist = dist; jetB->NN = jetA;}
2708 : }
2709 : }
2710 : }
2711 : }
2712 0 : vector<double> diJs(n);
2713 0 : for (int i = 0; i < n; i++) {
2714 0 : diJs[i] = _bj_diJ(&briefjets[i]);
2715 0 : briefjets[i].label_minheap_update_done();
2716 : }
2717 0 : MinHeap minheap(diJs);
2718 0 : vector<TiledJet *> jets_for_minheap;
2719 0 : jets_for_minheap.reserve(n);
2720 0 : int history_location = n-1;
2721 0 : while (n > 0) {
2722 0 : double diJ_min = minheap.minval() *_invR2;
2723 0 : jetA = head + minheap.minloc();
2724 0 : history_location++;
2725 0 : jetB = jetA->NN;
2726 0 : if (jetB != NULL) {
2727 0 : if (jetA < jetB) {std::swap(jetA,jetB);}
2728 0 : int nn; // new jet index
2729 0 : _do_ij_recombination_step(jetA->_jets_index, jetB->_jets_index, diJ_min, nn);
2730 0 : _bj_remove_from_tiles(jetA);
2731 0 : oldB = * jetB; // take a copy because we will need it...
2732 0 : _bj_remove_from_tiles(jetB);
2733 0 : _tj_set_jetinfo(jetB, nn); // cause jetB to become _jets[nn]
2734 0 : } else {
2735 0 : _do_iB_recombination_step(jetA->_jets_index, diJ_min);
2736 0 : _bj_remove_from_tiles(jetA);
2737 : }
2738 0 : minheap.remove(jetA-head);
2739 0 : int n_near_tiles = 0;
2740 0 : _add_untagged_neighbours_to_tile_union(jetA->tile_index,
2741 : tile_union, n_near_tiles);
2742 0 : if (jetB != NULL) {
2743 0 : if (jetB->tile_index != jetA->tile_index) {
2744 0 : _add_untagged_neighbours_to_tile_union(jetB->tile_index,
2745 : tile_union,n_near_tiles);
2746 : }
2747 0 : if (oldB.tile_index != jetA->tile_index &&
2748 0 : oldB.tile_index != jetB->tile_index) {
2749 : // GS: the line below generates a warning that oldB.tile_index
2750 : // may be used uninitialised. However, to reach this point, we
2751 : // ned jetB != NULL (see test a few lines above) and is jetB
2752 : // !=NULL, one would have gone through "oldB = *jetB before
2753 : // (see piece of code ~20 line above), so the index is
2754 : // initialised. We do not do anything to avoid the warning to
2755 : // avoid any potential speed impact.
2756 0 : _add_untagged_neighbours_to_tile_union(oldB.tile_index,
2757 : tile_union,n_near_tiles);
2758 : }
2759 0 : jetB->label_minheap_update_needed();
2760 0 : jets_for_minheap.push_back(jetB);
2761 : }
2762 0 : for (int itile = 0; itile < n_near_tiles; itile++) {
2763 0 : Tile * tile_ptr = &_tiles[tile_union[itile]];
2764 0 : tile_ptr->tagged = false; // reset tag, since we're done with unions
2765 0 : for (TiledJet * jetI = tile_ptr->head; jetI != NULL; jetI = jetI->next) {
2766 : // see if jetI had jetA or jetB as a NN -- if so recalculate the NN
2767 0 : if (jetI->NN == jetA || (jetI->NN == jetB && jetB != NULL)) {
2768 0 : jetI->NN_dist = _R2;
2769 0 : jetI->NN = NULL;
2770 : // label jetI as needing heap action...
2771 0 : if (!jetI->minheap_update_needed()) {
2772 0 : jetI->label_minheap_update_needed();
2773 0 : jets_for_minheap.push_back(jetI);}
2774 : // now go over tiles that are neighbours of I (include own tile)
2775 0 : for (Tile ** near_tile = tile_ptr->begin_tiles;
2776 0 : near_tile != tile_ptr->end_tiles; near_tile++) {
2777 : // and then over the contents of that tile
2778 0 : for (TiledJet * jetJ = (*near_tile)->head;
2779 0 : jetJ != NULL; jetJ = jetJ->next) {
2780 0 : double dist = _bj_dist(jetI,jetJ);
2781 0 : if (dist < jetI->NN_dist && jetJ != jetI) {
2782 0 : jetI->NN_dist = dist; jetI->NN = jetJ;
2783 0 : }
2784 : }
2785 : }
2786 0 : }
2787 : // check whether new jetB is closer than jetI's current NN and
2788 : // if jetI is closer than jetB's current (evolving) nearest
2789 : // neighbour. Where relevant update things
2790 0 : if (jetB != NULL) {
2791 0 : double dist = _bj_dist(jetI,jetB);
2792 0 : if (dist < jetI->NN_dist) {
2793 0 : if (jetI != jetB) {
2794 0 : jetI->NN_dist = dist;
2795 0 : jetI->NN = jetB;
2796 : // label jetI as needing heap action...
2797 0 : if (!jetI->minheap_update_needed()) {
2798 0 : jetI->label_minheap_update_needed();
2799 0 : jets_for_minheap.push_back(jetI);}
2800 : }
2801 : }
2802 0 : if (dist < jetB->NN_dist) {
2803 0 : if (jetI != jetB) {
2804 0 : jetB->NN_dist = dist;
2805 0 : jetB->NN = jetI;}
2806 : }
2807 0 : }
2808 : }
2809 : }
2810 0 : while (jets_for_minheap.size() > 0) {
2811 0 : TiledJet * jetI = jets_for_minheap.back();
2812 0 : jets_for_minheap.pop_back();
2813 0 : minheap.update(jetI-head, _bj_diJ(jetI));
2814 0 : jetI->label_minheap_update_done();
2815 : }
2816 0 : n--;
2817 0 : }
2818 0 : delete[] briefjets;
2819 0 : }
2820 : FJCORE_END_NAMESPACE
2821 : FJCORE_BEGIN_NAMESPACE // defined in fastjet/internal/base.hh
2822 : using namespace std;
2823 0 : CompositeJetStructure::CompositeJetStructure(const std::vector<PseudoJet> & initial_pieces,
2824 : const JetDefinition::Recombiner * recombiner)
2825 0 : : _pieces(initial_pieces){
2826 : if (recombiner){}; // ugly trick to prevent a gcc warning
2827 0 : _area_4vector_ptr = 0;
2828 0 : }
2829 : std::string CompositeJetStructure::description() const{
2830 0 : string str = "Composite PseudoJet";
2831 : return str;
2832 0 : }
2833 : bool CompositeJetStructure::has_constituents() const{
2834 0 : return _pieces.size()!=0;
2835 : }
2836 : std::vector<PseudoJet> CompositeJetStructure::constituents(const PseudoJet & /*jet*/) const{
2837 0 : vector<PseudoJet> all_constituents;
2838 0 : for (unsigned i = 0; i < _pieces.size(); i++) {
2839 0 : if (_pieces[i].has_constituents()){
2840 0 : vector<PseudoJet> constits = _pieces[i].constituents();
2841 0 : copy(constits.begin(), constits.end(), back_inserter(all_constituents));
2842 0 : } else {
2843 0 : all_constituents.push_back(_pieces[i]);
2844 : }
2845 : }
2846 : return all_constituents;
2847 0 : }
2848 : std::vector<PseudoJet> CompositeJetStructure::pieces(const PseudoJet & /*jet*/) const{
2849 0 : return _pieces;
2850 : }
2851 : FJCORE_END_NAMESPACE // defined in fastjet/internal/base.hh
2852 : #include <sstream>
2853 : #ifdef FJCORE_HAVE_EXECINFO_H
2854 : #include <execinfo.h>
2855 : #include <cstdlib>
2856 : #endif
2857 : FJCORE_BEGIN_NAMESPACE // defined in fastjet/internal/base.hh
2858 : using namespace std;
2859 : bool Error::_print_errors = true;
2860 : bool Error::_print_backtrace = false;
2861 : ostream * Error::_default_ostr = & cerr;
2862 0 : Error::Error(const std::string & message_in) {
2863 0 : _message = message_in;
2864 0 : if (_print_errors && _default_ostr){
2865 0 : ostringstream oss;
2866 0 : oss << "fjcore::Error: "<< message_in << endl;
2867 : #ifdef FJCORE_HAVE_EXECINFO_H
2868 0 : if (_print_backtrace){
2869 0 : void * array[10];
2870 : char ** messages;
2871 0 : int size = backtrace(array, 10);
2872 0 : messages = backtrace_symbols(array, size);
2873 0 : oss << "stack:" << endl;
2874 0 : for (int i = 1; i < size && messages != NULL; ++i){
2875 0 : oss << " #" << i << ": " << messages[i] << endl;
2876 : }
2877 0 : free(messages);
2878 0 : }
2879 : #endif
2880 0 : *_default_ostr << oss.str();
2881 0 : _default_ostr->flush();
2882 0 : }
2883 0 : }
2884 : FJCORE_END_NAMESPACE
2885 : #include <string>
2886 : #include <sstream>
2887 : using namespace std;
2888 : FJCORE_BEGIN_NAMESPACE
2889 : FJCORE_END_NAMESPACE
2890 : #include<sstream>
2891 : FJCORE_BEGIN_NAMESPACE // defined in fastjet/internal/base.hh
2892 : using namespace std;
2893 : const double JetDefinition::max_allowable_R = 1000.0;
2894 0 : JetDefinition::JetDefinition(JetAlgorithm jet_algorithm_in,
2895 : double R_in,
2896 : Strategy strategy_in,
2897 : RecombinationScheme recomb_scheme_in,
2898 : int nparameters) :
2899 0 : _jet_algorithm(jet_algorithm_in), _Rparam(R_in), _strategy(strategy_in) {
2900 0 : if (_jet_algorithm == ee_kt_algorithm) {
2901 0 : _Rparam = 4.0; // introduce a fictional R that ensures that
2902 0 : } else {
2903 0 : if (R_in > max_allowable_R) {
2904 0 : ostringstream oss;
2905 0 : oss << "Requested R = " << R_in << " for jet definition is larger than max_allowable_R = " << max_allowable_R;
2906 0 : throw Error(oss.str());
2907 0 : }
2908 : }
2909 0 : switch (jet_algorithm_in) {
2910 : case ee_kt_algorithm:
2911 0 : if (nparameters != 0) {
2912 0 : ostringstream oss;
2913 0 : oss << "ee_kt_algorithm should be constructed with 0 parameters but was called with "
2914 0 : << nparameters << " parameter(s)\n";
2915 0 : throw Error(oss.str());
2916 0 : }
2917 : break;
2918 : case genkt_algorithm:
2919 : case ee_genkt_algorithm:
2920 0 : if (nparameters != 2) {
2921 0 : ostringstream oss;
2922 0 : oss << "(ee_)genkt_algorithm should be constructed with 2 parameters but was called with "
2923 0 : << nparameters << " parameter(s)\n";
2924 0 : throw Error(oss.str());
2925 0 : }
2926 : break;
2927 : default:
2928 0 : if (nparameters != 1) {
2929 0 : ostringstream oss;
2930 0 : oss << "The jet algorithm you requested ("
2931 0 : << jet_algorithm_in << ") should be constructed with 1 parameter but was called with "
2932 0 : << nparameters << " parameter(s)\n";
2933 0 : throw Error(oss.str());
2934 0 : }
2935 : }
2936 0 : assert (_strategy != plugin_strategy);
2937 0 : _plugin = NULL;
2938 0 : set_recombination_scheme(recomb_scheme_in);
2939 0 : set_extra_param(0.0); // make sure it's defined
2940 0 : }
2941 : string JetDefinition::description() const {
2942 0 : ostringstream name;
2943 0 : if (jet_algorithm() == plugin_algorithm) {
2944 0 : return plugin()->description();
2945 0 : } else if (jet_algorithm() == kt_algorithm) {
2946 0 : name << "Longitudinally invariant kt algorithm with R = " << R();
2947 0 : name << " and " << recombiner()->description();
2948 0 : } else if (jet_algorithm() == cambridge_algorithm) {
2949 0 : name << "Longitudinally invariant Cambridge/Aachen algorithm with R = "
2950 0 : << R() ;
2951 0 : name << " and " << recombiner()->description();
2952 0 : } else if (jet_algorithm() == antikt_algorithm) {
2953 0 : name << "Longitudinally invariant anti-kt algorithm with R = "
2954 0 : << R() ;
2955 0 : name << " and " << recombiner()->description();
2956 0 : } else if (jet_algorithm() == genkt_algorithm) {
2957 0 : name << "Longitudinally invariant generalised kt algorithm with R = "
2958 0 : << R() << ", p = " << extra_param();
2959 0 : name << " and " << recombiner()->description();
2960 0 : } else if (jet_algorithm() == cambridge_for_passive_algorithm) {
2961 0 : name << "Longitudinally invariant Cambridge/Aachen algorithm with R = "
2962 0 : << R() << "and a special hack whereby particles with kt < "
2963 0 : << extra_param() << "are treated as passive ghosts";
2964 0 : } else if (jet_algorithm() == ee_kt_algorithm) {
2965 0 : name << "e+e- kt (Durham) algorithm (NB: no R)";
2966 0 : name << " with " << recombiner()->description();
2967 0 : } else if (jet_algorithm() == ee_genkt_algorithm) {
2968 0 : name << "e+e- generalised kt algorithm with R = "
2969 0 : << R() << ", p = " << extra_param();
2970 0 : name << " and " << recombiner()->description();
2971 0 : } else if (jet_algorithm() == undefined_jet_algorithm) {
2972 0 : name << "uninitialised JetDefinition (jet_algorithm=undefined_jet_algorithm)" ;
2973 : } else {
2974 0 : throw Error("JetDefinition::description(): unrecognized jet_algorithm");
2975 : }
2976 0 : return name.str();
2977 0 : }
2978 : void JetDefinition::set_recombination_scheme(
2979 : RecombinationScheme recomb_scheme) {
2980 0 : _default_recombiner = JetDefinition::DefaultRecombiner(recomb_scheme);
2981 0 : if (_recombiner_shared()) _recombiner_shared.reset();
2982 0 : _recombiner = 0;
2983 0 : }
2984 : bool JetDefinition::has_same_recombiner(const JetDefinition &other_jd) const{
2985 0 : const RecombinationScheme & scheme = recombination_scheme();
2986 0 : if (other_jd.recombination_scheme() != scheme) return false;
2987 0 : return (scheme != external_scheme)
2988 0 : || (recombiner() == other_jd.recombiner());
2989 0 : }
2990 : void JetDefinition::delete_recombiner_when_unused(){
2991 0 : if (_recombiner == 0){
2992 0 : throw Error("tried to call JetDefinition::delete_recombiner_when_unused() for a JetDefinition without a user-defined recombination scheme");
2993 : }
2994 0 : _recombiner_shared.reset(_recombiner);
2995 0 : }
2996 : void JetDefinition::delete_plugin_when_unused(){
2997 0 : if (_plugin == 0){
2998 0 : throw Error("tried to call JetDefinition::delete_plugin_when_unused() for a JetDefinition without a plugin");
2999 : }
3000 0 : _plugin_shared.reset(_plugin);
3001 0 : }
3002 : string JetDefinition::DefaultRecombiner::description() const {
3003 0 : switch(_recomb_scheme) {
3004 : case E_scheme:
3005 0 : return "E scheme recombination";
3006 : case pt_scheme:
3007 0 : return "pt scheme recombination";
3008 : case pt2_scheme:
3009 0 : return "pt2 scheme recombination";
3010 : case Et_scheme:
3011 0 : return "Et scheme recombination";
3012 : case Et2_scheme:
3013 0 : return "Et2 scheme recombination";
3014 : case BIpt_scheme:
3015 0 : return "boost-invariant pt scheme recombination";
3016 : case BIpt2_scheme:
3017 0 : return "boost-invariant pt2 scheme recombination";
3018 : default:
3019 0 : ostringstream err;
3020 0 : err << "DefaultRecombiner: unrecognized recombination scheme "
3021 0 : << _recomb_scheme;
3022 0 : throw Error(err.str());
3023 0 : }
3024 0 : }
3025 : void JetDefinition::DefaultRecombiner::recombine(
3026 : const PseudoJet & pa, const PseudoJet & pb,
3027 : PseudoJet & pab) const {
3028 : double weighta, weightb;
3029 0 : switch(_recomb_scheme) {
3030 : case E_scheme:
3031 0 : pab.reset(pa.px()+pb.px(),
3032 0 : pa.py()+pb.py(),
3033 0 : pa.pz()+pb.pz(),
3034 0 : pa.E ()+pb.E ());
3035 0 : return;
3036 : case pt_scheme:
3037 : case Et_scheme:
3038 : case BIpt_scheme:
3039 0 : weighta = pa.perp();
3040 0 : weightb = pb.perp();
3041 0 : break;
3042 : case pt2_scheme:
3043 : case Et2_scheme:
3044 : case BIpt2_scheme:
3045 0 : weighta = pa.perp2();
3046 0 : weightb = pb.perp2();
3047 0 : break;
3048 : default:
3049 0 : ostringstream err;
3050 0 : err << "DefaultRecombiner: unrecognized recombination scheme "
3051 0 : << _recomb_scheme;
3052 0 : throw Error(err.str());
3053 0 : }
3054 0 : double perp_ab = pa.perp() + pb.perp();
3055 0 : if (perp_ab != 0.0) { // weights also non-zero...
3056 0 : double y_ab = (weighta * pa.rap() + weightb * pb.rap())/(weighta+weightb);
3057 0 : double phi_a = pa.phi(), phi_b = pb.phi();
3058 0 : if (phi_a - phi_b > pi) phi_b += twopi;
3059 0 : if (phi_a - phi_b < -pi) phi_b -= twopi;
3060 0 : double phi_ab = (weighta * phi_a + weightb * phi_b)/(weighta+weightb);
3061 0 : pab.reset_PtYPhiM(perp_ab,y_ab,phi_ab);
3062 0 : } else { // weights are zero
3063 0 : pab.reset(0.0, 0.0, 0.0, 0.0);
3064 : }
3065 0 : }
3066 : void JetDefinition::DefaultRecombiner::preprocess(PseudoJet & p) const {
3067 0 : switch(_recomb_scheme) {
3068 : case E_scheme:
3069 : case BIpt_scheme:
3070 : case BIpt2_scheme:
3071 : break;
3072 : case pt_scheme:
3073 : case pt2_scheme:
3074 : {
3075 0 : double newE = sqrt(p.perp2()+p.pz()*p.pz());
3076 0 : p.reset_momentum(p.px(), p.py(), p.pz(), newE);
3077 : }
3078 0 : break;
3079 : case Et_scheme:
3080 : case Et2_scheme:
3081 : {
3082 0 : double rescale = p.E()/sqrt(p.perp2()+p.pz()*p.pz());
3083 0 : p.reset_momentum(rescale*p.px(), rescale*p.py(), rescale*p.pz(), p.E());
3084 : }
3085 0 : break;
3086 : default:
3087 0 : ostringstream err;
3088 0 : err << "DefaultRecombiner: unrecognized recombination scheme "
3089 0 : << _recomb_scheme;
3090 0 : throw Error(err.str());
3091 0 : }
3092 0 : }
3093 : void JetDefinition::Plugin::set_ghost_separation_scale(double /*scale*/) const {
3094 0 : throw Error("set_ghost_separation_scale not supported");
3095 0 : }
3096 : PseudoJet join(const vector<PseudoJet> & pieces, const JetDefinition::Recombiner & recombiner){
3097 0 : PseudoJet result; // automatically initialised to 0
3098 0 : if (pieces.size()>0){
3099 0 : result = pieces[0];
3100 0 : for (unsigned int i=1; i<pieces.size(); i++)
3101 0 : recombiner.plus_equal(result, pieces[i]);
3102 0 : }
3103 0 : CompositeJetStructure *cj_struct = new CompositeJetStructure(pieces, &recombiner);
3104 0 : result.set_structure_shared_ptr(SharedPtr<PseudoJetStructureBase>(cj_struct));
3105 : return result;
3106 0 : }
3107 : PseudoJet join(const PseudoJet & j1,
3108 : const JetDefinition::Recombiner & recombiner){
3109 0 : return join(vector<PseudoJet>(1,j1), recombiner);
3110 0 : }
3111 : PseudoJet join(const PseudoJet & j1, const PseudoJet & j2,
3112 : const JetDefinition::Recombiner & recombiner){
3113 0 : vector<PseudoJet> pieces;
3114 0 : pieces.push_back(j1);
3115 0 : pieces.push_back(j2);
3116 0 : return join(pieces, recombiner);
3117 0 : }
3118 : PseudoJet join(const PseudoJet & j1, const PseudoJet & j2, const PseudoJet & j3,
3119 : const JetDefinition::Recombiner & recombiner){
3120 0 : vector<PseudoJet> pieces;
3121 0 : pieces.push_back(j1);
3122 0 : pieces.push_back(j2);
3123 0 : pieces.push_back(j3);
3124 0 : return join(pieces, recombiner);
3125 0 : }
3126 : PseudoJet join(const PseudoJet & j1, const PseudoJet & j2, const PseudoJet & j3, const PseudoJet & j4,
3127 : const JetDefinition::Recombiner & recombiner){
3128 0 : vector<PseudoJet> pieces;
3129 0 : pieces.push_back(j1);
3130 0 : pieces.push_back(j2);
3131 0 : pieces.push_back(j3);
3132 0 : pieces.push_back(j4);
3133 0 : return join(pieces, recombiner);
3134 0 : }
3135 : FJCORE_END_NAMESPACE
3136 : #include <sstream>
3137 : #include <limits>
3138 : using namespace std;
3139 : FJCORE_BEGIN_NAMESPACE
3140 : ostream * LimitedWarning::_default_ostr = &cerr;
3141 6 : std::list< LimitedWarning::Summary > LimitedWarning::_global_warnings_summary;
3142 : int LimitedWarning::_max_warn_default = 5;
3143 : void LimitedWarning::warn(const std::string & warning) {
3144 0 : warn(warning, _default_ostr);
3145 0 : }
3146 : void LimitedWarning::warn(const std::string & warning, std::ostream * ostr) {
3147 0 : if (_this_warning_summary == 0) {
3148 0 : _global_warnings_summary.push_back(Summary(warning, 0));
3149 0 : _this_warning_summary = & (_global_warnings_summary.back());
3150 0 : }
3151 0 : if (_n_warn_so_far < _max_warn) {
3152 0 : ostringstream warnstr;
3153 0 : warnstr << "WARNING: ";
3154 0 : warnstr << warning;
3155 0 : _n_warn_so_far++;
3156 0 : if (_n_warn_so_far == _max_warn) warnstr << " (LAST SUCH WARNING)";
3157 0 : warnstr << std::endl;
3158 0 : if (ostr) {
3159 0 : (*ostr) << warnstr.str();
3160 0 : ostr->flush(); // get something written to file even if the program aborts
3161 : }
3162 0 : }
3163 0 : if (_this_warning_summary->second < numeric_limits<unsigned>::max()) {
3164 0 : _this_warning_summary->second++;
3165 0 : }
3166 0 : }
3167 : string LimitedWarning::summary() {
3168 0 : ostringstream str;
3169 0 : for (list<Summary>::const_iterator it = _global_warnings_summary.begin();
3170 0 : it != _global_warnings_summary.end(); it++) {
3171 0 : str << it->second << " times: " << it->first << endl;
3172 : }
3173 0 : return str.str();
3174 0 : }
3175 : FJCORE_END_NAMESPACE
3176 : #include<iostream>
3177 : #include<cmath>
3178 : #include<limits>
3179 : FJCORE_BEGIN_NAMESPACE // defined in fastjet/internal/base.hh
3180 : using namespace std;
3181 : void MinHeap::_initialise(const std::vector<double> & values){
3182 0 : for (unsigned i = values.size(); i < _heap.size(); i++) {
3183 0 : _heap[i].value = std::numeric_limits<double>::max();
3184 0 : _heap[i].minloc = &(_heap[i]);
3185 : }
3186 0 : for (unsigned i = 0; i < values.size(); i++) {
3187 0 : _heap[i].value = values[i];
3188 0 : _heap[i].minloc = &(_heap[i]);
3189 : }
3190 0 : for (unsigned i = _heap.size()-1; i > 0; i--) {
3191 0 : ValueLoc * parent = &(_heap[(i-1)/2]);
3192 0 : ValueLoc * here = &(_heap[i]);
3193 0 : if (here->minloc->value < parent->minloc->value) {
3194 0 : parent->minloc = here->minloc;
3195 0 : }
3196 : }
3197 0 : }
3198 : void MinHeap::update(unsigned int loc, double new_value) {
3199 0 : assert(loc < _heap.size());
3200 0 : ValueLoc * start = &(_heap[loc]);
3201 0 : if (start->minloc != start && !(new_value < start->minloc->value)) {
3202 0 : start->value = new_value;
3203 0 : return;
3204 : }
3205 0 : start->value = new_value;
3206 0 : start->minloc = start;
3207 : bool change_made = true;
3208 0 : ValueLoc * heap_end = (&(_heap[0])) + _heap.size();
3209 0 : while(change_made) {
3210 0 : ValueLoc * here = &(_heap[loc]);
3211 : change_made = false;
3212 0 : if (here->minloc == start) {
3213 0 : here->minloc = here; change_made = true;
3214 0 : }
3215 0 : ValueLoc * child = &(_heap[2*loc+1]);
3216 0 : if (child < heap_end && child->minloc->value < here->minloc->value ) {
3217 0 : here->minloc = child->minloc;
3218 0 : change_made = true;}
3219 0 : child++;
3220 0 : if (child < heap_end && child->minloc->value < here->minloc->value ) {
3221 0 : here->minloc = child->minloc;
3222 0 : change_made = true;}
3223 0 : if (loc == 0) {break;}
3224 0 : loc = (loc-1)/2;
3225 0 : }
3226 0 : }
3227 : FJCORE_END_NAMESPACE
3228 : #include<valarray>
3229 : #include<iostream>
3230 : #include<sstream>
3231 : #include<cmath>
3232 : #include<algorithm>
3233 : #include <cstdarg>
3234 : FJCORE_BEGIN_NAMESPACE // defined in fastjet/internal/base.hh
3235 : using namespace std;
3236 0 : PseudoJet::PseudoJet(const double px_in, const double py_in, const double pz_in, const double E_in) {
3237 0 : _E = E_in ;
3238 0 : _px = px_in;
3239 0 : _py = py_in;
3240 0 : _pz = pz_in;
3241 0 : this->_finish_init();
3242 0 : _reset_indices();
3243 0 : }
3244 : void PseudoJet::_finish_init () {
3245 0 : _kt2 = this->px()*this->px() + this->py()*this->py();
3246 0 : _phi = pseudojet_invalid_phi;
3247 0 : _rap = pseudojet_invalid_rap;
3248 0 : }
3249 : void PseudoJet::_set_rap_phi() const {
3250 0 : if (_kt2 == 0.0) {
3251 0 : _phi = 0.0; }
3252 : else {
3253 0 : _phi = atan2(this->py(),this->px());
3254 : }
3255 0 : if (_phi < 0.0) {_phi += twopi;}
3256 0 : if (_phi >= twopi) {_phi -= twopi;} // can happen if phi=-|eps<1e-15|?
3257 0 : if (this->E() == abs(this->pz()) && _kt2 == 0) {
3258 0 : double MaxRapHere = MaxRap + abs(this->pz());
3259 0 : if (this->pz() >= 0.0) {_rap = MaxRapHere;} else {_rap = -MaxRapHere;}
3260 0 : } else {
3261 0 : double effective_m2 = max(0.0,m2()); // force non tachyonic mass
3262 0 : double E_plus_pz = _E + abs(_pz); // the safer of p+, p-
3263 0 : _rap = 0.5*log((_kt2 + effective_m2)/(E_plus_pz*E_plus_pz));
3264 0 : if (_pz > 0) {_rap = - _rap;}
3265 : }
3266 0 : }
3267 : valarray<double> PseudoJet::four_mom() const {
3268 0 : valarray<double> mom(4);
3269 0 : mom[0] = _px;
3270 0 : mom[1] = _py;
3271 0 : mom[2] = _pz;
3272 0 : mom[3] = _E ;
3273 : return mom;
3274 0 : }
3275 : double PseudoJet::operator () (int i) const {
3276 0 : switch(i) {
3277 : case X:
3278 0 : return px();
3279 : case Y:
3280 0 : return py();
3281 : case Z:
3282 0 : return pz();
3283 : case T:
3284 0 : return e();
3285 : default:
3286 0 : ostringstream err;
3287 0 : err << "PseudoJet subscripting: bad index (" << i << ")";
3288 0 : throw Error(err.str());
3289 0 : }
3290 : return 0.;
3291 0 : }
3292 : double PseudoJet::pseudorapidity() const {
3293 0 : if (px() == 0.0 && py() ==0.0) return MaxRap;
3294 0 : if (pz() == 0.0) return 0.0;
3295 0 : double theta = atan(perp()/pz());
3296 0 : if (theta < 0) theta += pi;
3297 0 : return -log(tan(theta/2));
3298 0 : }
3299 : PseudoJet operator+ (const PseudoJet & jet1, const PseudoJet & jet2) {
3300 0 : return PseudoJet(jet1.px()+jet2.px(),
3301 0 : jet1.py()+jet2.py(),
3302 0 : jet1.pz()+jet2.pz(),
3303 0 : jet1.E() +jet2.E() );
3304 : }
3305 : PseudoJet operator- (const PseudoJet & jet1, const PseudoJet & jet2) {
3306 0 : return PseudoJet(jet1.px()-jet2.px(),
3307 0 : jet1.py()-jet2.py(),
3308 0 : jet1.pz()-jet2.pz(),
3309 0 : jet1.E() -jet2.E() );
3310 : }
3311 : PseudoJet operator* (double coeff, const PseudoJet & jet) {
3312 0 : PseudoJet coeff_times_jet(jet);
3313 0 : coeff_times_jet *= coeff;
3314 : return coeff_times_jet;
3315 0 : }
3316 : PseudoJet operator* (const PseudoJet & jet, double coeff) {
3317 0 : return coeff*jet;
3318 : }
3319 : PseudoJet operator/ (const PseudoJet & jet, double coeff) {
3320 0 : return (1.0/coeff)*jet;
3321 : }
3322 : void PseudoJet::operator*=(double coeff) {
3323 0 : _px *= coeff;
3324 0 : _py *= coeff;
3325 0 : _pz *= coeff;
3326 0 : _E *= coeff;
3327 0 : _kt2*= coeff*coeff;
3328 0 : }
3329 : void PseudoJet::operator/=(double coeff) {
3330 0 : (*this) *= 1.0/coeff;
3331 0 : }
3332 : void PseudoJet::operator+=(const PseudoJet & other_jet) {
3333 0 : _px += other_jet._px;
3334 0 : _py += other_jet._py;
3335 0 : _pz += other_jet._pz;
3336 0 : _E += other_jet._E ;
3337 0 : _finish_init(); // we need to recalculate phi,rap,kt2
3338 0 : }
3339 : void PseudoJet::operator-=(const PseudoJet & other_jet) {
3340 0 : _px -= other_jet._px;
3341 0 : _py -= other_jet._py;
3342 0 : _pz -= other_jet._pz;
3343 0 : _E -= other_jet._E ;
3344 0 : _finish_init(); // we need to recalculate phi,rap,kt2
3345 0 : }
3346 : bool operator==(const PseudoJet & a, const PseudoJet & b) {
3347 0 : if (a.px() != b.px()) return false;
3348 0 : if (a.py() != b.py()) return false;
3349 0 : if (a.pz() != b.pz()) return false;
3350 0 : if (a.E () != b.E ()) return false;
3351 0 : if (a.user_index() != b.user_index()) return false;
3352 0 : if (a.cluster_hist_index() != b.cluster_hist_index()) return false;
3353 0 : if (a.user_info_ptr() != b.user_info_ptr()) return false;
3354 0 : if (a.structure_ptr() != b.structure_ptr()) return false;
3355 0 : return true;
3356 0 : }
3357 : bool operator==(const PseudoJet & jet, const double val) {
3358 0 : if (val != 0)
3359 0 : throw Error("comparing a PseudoJet with a non-zero constant (double) is not allowed.");
3360 0 : return (jet.px() == 0 && jet.py() == 0 &&
3361 0 : jet.pz() == 0 && jet.E() == 0);
3362 0 : }
3363 : PseudoJet & PseudoJet::boost(const PseudoJet & prest) {
3364 0 : if (prest.px() == 0.0 && prest.py() == 0.0 && prest.pz() == 0.0)
3365 0 : return *this;
3366 0 : double m_local = prest.m();
3367 0 : assert(m_local != 0);
3368 0 : double pf4 = ( px()*prest.px() + py()*prest.py()
3369 0 : + pz()*prest.pz() + E()*prest.E() )/m_local;
3370 0 : double fn = (pf4 + E()) / (prest.E() + m_local);
3371 0 : _px += fn*prest.px();
3372 0 : _py += fn*prest.py();
3373 0 : _pz += fn*prest.pz();
3374 0 : _E = pf4;
3375 0 : _finish_init(); // we need to recalculate phi,rap,kt2
3376 : return *this;
3377 0 : }
3378 : PseudoJet & PseudoJet::unboost(const PseudoJet & prest) {
3379 0 : if (prest.px() == 0.0 && prest.py() == 0.0 && prest.pz() == 0.0)
3380 0 : return *this;
3381 0 : double m_local = prest.m();
3382 0 : assert(m_local != 0);
3383 0 : double pf4 = ( -px()*prest.px() - py()*prest.py()
3384 0 : - pz()*prest.pz() + E()*prest.E() )/m_local;
3385 0 : double fn = (pf4 + E()) / (prest.E() + m_local);
3386 0 : _px -= fn*prest.px();
3387 0 : _py -= fn*prest.py();
3388 0 : _pz -= fn*prest.pz();
3389 0 : _E = pf4;
3390 0 : _finish_init(); // we need to recalculate phi,rap,kt2
3391 : return *this;
3392 0 : }
3393 : bool have_same_momentum(const PseudoJet & jeta, const PseudoJet & jetb) {
3394 0 : return jeta.px() == jetb.px()
3395 0 : && jeta.py() == jetb.py()
3396 0 : && jeta.pz() == jetb.pz()
3397 0 : && jeta.E() == jetb.E();
3398 : }
3399 : void PseudoJet::set_cached_rap_phi(double rap_in, double phi_in) {
3400 0 : _rap = rap_in; _phi = phi_in;
3401 0 : if (_phi >= twopi) _phi -= twopi;
3402 0 : if (_phi < 0) _phi += twopi;
3403 0 : }
3404 : void PseudoJet::reset_momentum_PtYPhiM(double pt_in, double y_in, double phi_in, double m_in) {
3405 0 : assert(phi_in < 2*twopi && phi_in > -twopi);
3406 0 : double ptm = (m_in == 0) ? pt_in : sqrt(pt_in*pt_in+m_in*m_in);
3407 0 : double exprap = exp(y_in);
3408 0 : double pminus = ptm/exprap;
3409 0 : double pplus = ptm*exprap;
3410 0 : double px_local = pt_in*cos(phi_in);
3411 0 : double py_local = pt_in*sin(phi_in);
3412 0 : reset_momentum(px_local,py_local,0.5*(pplus-pminus),0.5*(pplus+pminus));
3413 0 : set_cached_rap_phi(y_in,phi_in);
3414 0 : }
3415 : PseudoJet PtYPhiM(double pt, double y, double phi, double m) {
3416 0 : assert(phi < 2*twopi && phi > -twopi);
3417 0 : double ptm = (m == 0) ? pt : sqrt(pt*pt+m*m);
3418 0 : double exprap = exp(y);
3419 0 : double pminus = ptm/exprap;
3420 0 : double pplus = ptm*exprap;
3421 0 : double px = pt*cos(phi);
3422 0 : double py = pt*sin(phi);
3423 0 : PseudoJet mom(px,py,0.5*(pplus-pminus),0.5*(pplus+pminus));
3424 0 : mom.set_cached_rap_phi(y,phi);
3425 : return mom;
3426 0 : }
3427 : double PseudoJet::kt_distance(const PseudoJet & other) const {
3428 0 : double distance = min(_kt2, other._kt2);
3429 0 : double dphi = abs(phi() - other.phi());
3430 0 : if (dphi > pi) {dphi = twopi - dphi;}
3431 0 : double drap = rap() - other.rap();
3432 0 : distance = distance * (dphi*dphi + drap*drap);
3433 0 : return distance;
3434 : }
3435 : double PseudoJet::plain_distance(const PseudoJet & other) const {
3436 0 : double dphi = abs(phi() - other.phi());
3437 0 : if (dphi > pi) {dphi = twopi - dphi;}
3438 0 : double drap = rap() - other.rap();
3439 0 : return (dphi*dphi + drap*drap);
3440 : }
3441 : double PseudoJet::delta_phi_to(const PseudoJet & other) const {
3442 0 : double dphi = other.phi() - phi();
3443 0 : if (dphi > pi) dphi -= twopi;
3444 0 : if (dphi < -pi) dphi += twopi;
3445 0 : return dphi;
3446 : }
3447 : string PseudoJet::description() const{
3448 0 : if (!_structure())
3449 0 : return "standard PseudoJet (with no associated clustering information)";
3450 0 : return _structure()->description();
3451 0 : }
3452 : bool PseudoJet::has_associated_cluster_sequence() const{
3453 0 : return (_structure()) && (_structure->has_associated_cluster_sequence());
3454 : }
3455 : const ClusterSequence* PseudoJet::associated_cluster_sequence() const{
3456 0 : if (! has_associated_cluster_sequence()) return NULL;
3457 0 : return _structure->associated_cluster_sequence();
3458 0 : }
3459 : bool PseudoJet::has_valid_cluster_sequence() const{
3460 0 : return (_structure()) && (_structure->has_valid_cluster_sequence());
3461 : }
3462 : const ClusterSequence * PseudoJet::validated_cs() const {
3463 0 : return validated_structure_ptr()->validated_cs();
3464 : }
3465 : void PseudoJet::set_structure_shared_ptr(const SharedPtr<PseudoJetStructureBase> &structure_in){
3466 0 : _structure = structure_in;
3467 0 : }
3468 : bool PseudoJet::has_structure() const{
3469 0 : return _structure();
3470 : }
3471 : const PseudoJetStructureBase* PseudoJet::structure_ptr() const {
3472 0 : if (!_structure()) return NULL;
3473 0 : return _structure();
3474 0 : }
3475 : PseudoJetStructureBase* PseudoJet::structure_non_const_ptr(){
3476 0 : if (!_structure()) return NULL;
3477 0 : return _structure();
3478 0 : }
3479 : const PseudoJetStructureBase* PseudoJet::validated_structure_ptr() const {
3480 0 : if (!_structure())
3481 0 : throw Error("Trying to access the structure of a PseudoJet which has no associated structure");
3482 0 : return _structure();
3483 0 : }
3484 : const SharedPtr<PseudoJetStructureBase> & PseudoJet::structure_shared_ptr() const {
3485 0 : return _structure;
3486 : }
3487 : bool PseudoJet::has_partner(PseudoJet &partner) const{
3488 0 : return validated_structure_ptr()->has_partner(*this, partner);
3489 : }
3490 : bool PseudoJet::has_child(PseudoJet &child) const{
3491 0 : return validated_structure_ptr()->has_child(*this, child);
3492 : }
3493 : bool PseudoJet::has_parents(PseudoJet &parent1, PseudoJet &parent2) const{
3494 0 : return validated_structure_ptr()->has_parents(*this, parent1, parent2);
3495 : }
3496 : bool PseudoJet::contains(const PseudoJet &constituent) const{
3497 0 : return validated_structure_ptr()->object_in_jet(constituent, *this);
3498 : }
3499 : bool PseudoJet::is_inside(const PseudoJet &jet) const{
3500 0 : return validated_structure_ptr()->object_in_jet(*this, jet);
3501 : }
3502 : bool PseudoJet::has_constituents() const{
3503 0 : return (_structure()) && (_structure->has_constituents());
3504 : }
3505 : vector<PseudoJet> PseudoJet::constituents() const{
3506 0 : return validated_structure_ptr()->constituents(*this);
3507 : }
3508 : bool PseudoJet::has_exclusive_subjets() const{
3509 0 : return (_structure()) && (_structure->has_exclusive_subjets());
3510 : }
3511 : std::vector<PseudoJet> PseudoJet::exclusive_subjets (const double & dcut) const {
3512 0 : return validated_structure_ptr()->exclusive_subjets(*this, dcut);
3513 : }
3514 : int PseudoJet::n_exclusive_subjets(const double & dcut) const {
3515 0 : return validated_structure_ptr()->n_exclusive_subjets(*this, dcut);
3516 : }
3517 : std::vector<PseudoJet> PseudoJet::exclusive_subjets_up_to (int nsub) const {
3518 0 : return validated_structure_ptr()->exclusive_subjets_up_to(*this, nsub);
3519 : }
3520 : std::vector<PseudoJet> PseudoJet::exclusive_subjets (int nsub) const {
3521 0 : vector<PseudoJet> subjets = exclusive_subjets_up_to(nsub);
3522 0 : if (int(subjets.size()) < nsub) {
3523 0 : ostringstream err;
3524 0 : err << "Requested " << nsub << " exclusive subjets, but there were only "
3525 0 : << subjets.size() << " particles in the jet";
3526 0 : throw Error(err.str());
3527 0 : }
3528 : return subjets;
3529 0 : }
3530 : double PseudoJet::exclusive_subdmerge(int nsub) const {
3531 0 : return validated_structure_ptr()->exclusive_subdmerge(*this, nsub);
3532 : }
3533 : double PseudoJet::exclusive_subdmerge_max(int nsub) const {
3534 0 : return validated_structure_ptr()->exclusive_subdmerge_max(*this, nsub);
3535 : }
3536 : bool PseudoJet::has_pieces() const{
3537 0 : return ((_structure()) && (_structure->has_pieces(*this)));
3538 : }
3539 : std::vector<PseudoJet> PseudoJet::pieces() const{
3540 0 : return validated_structure_ptr()->pieces(*this);
3541 : }
3542 0 : PseudoJet::InexistentUserInfo::InexistentUserInfo() : Error("you attempted to perform a dynamic cast of a PseudoJet's extra info, but the extra info pointer was null")
3543 0 : {}
3544 : void sort_indices(vector<int> & indices,
3545 : const vector<double> & values) {
3546 0 : IndexedSortHelper index_sort_helper(&values);
3547 0 : sort(indices.begin(), indices.end(), index_sort_helper);
3548 0 : }
3549 : template<class T> vector<T> objects_sorted_by_values(
3550 : const vector<T> & objects,
3551 : const vector<double> & values) {
3552 0 : assert(objects.size() == values.size());
3553 0 : vector<int> indices(values.size());
3554 0 : for (size_t i = 0; i < indices.size(); i++) {indices[i] = i;}
3555 0 : sort_indices(indices, values);
3556 0 : vector<T> objects_sorted(objects.size());
3557 0 : for (size_t i = 0; i < indices.size(); i++) {
3558 0 : objects_sorted[i] = objects[indices[i]];
3559 : }
3560 : return objects_sorted;
3561 0 : }
3562 : vector<PseudoJet> sorted_by_pt(const vector<PseudoJet> & jets) {
3563 0 : vector<double> minus_kt2(jets.size());
3564 0 : for (size_t i = 0; i < jets.size(); i++) {minus_kt2[i] = -jets[i].kt2();}
3565 0 : return objects_sorted_by_values(jets, minus_kt2);
3566 0 : }
3567 : vector<PseudoJet> sorted_by_rapidity(const vector<PseudoJet> & jets) {
3568 0 : vector<double> rapidities(jets.size());
3569 0 : for (size_t i = 0; i < jets.size(); i++) {rapidities[i] = jets[i].rap();}
3570 0 : return objects_sorted_by_values(jets, rapidities);
3571 0 : }
3572 : vector<PseudoJet> sorted_by_E(const vector<PseudoJet> & jets) {
3573 0 : vector<double> energies(jets.size());
3574 0 : for (size_t i = 0; i < jets.size(); i++) {energies[i] = -jets[i].E();}
3575 0 : return objects_sorted_by_values(jets, energies);
3576 0 : }
3577 : vector<PseudoJet> sorted_by_pz(const vector<PseudoJet> & jets) {
3578 0 : vector<double> pz(jets.size());
3579 0 : for (size_t i = 0; i < jets.size(); i++) {pz[i] = jets[i].pz();}
3580 0 : return objects_sorted_by_values(jets, pz);
3581 0 : }
3582 : PseudoJet join(const vector<PseudoJet> & pieces){
3583 0 : PseudoJet result; // automatically initialised to 0
3584 0 : for (unsigned int i=0; i<pieces.size(); i++)
3585 0 : result += pieces[i];
3586 0 : CompositeJetStructure *cj_struct = new CompositeJetStructure(pieces);
3587 0 : result.set_structure_shared_ptr(SharedPtr<PseudoJetStructureBase>(cj_struct));
3588 : return result;
3589 0 : }
3590 : PseudoJet join(const PseudoJet & j1){
3591 0 : return join(vector<PseudoJet>(1,j1));
3592 0 : }
3593 : PseudoJet join(const PseudoJet & j1, const PseudoJet & j2){
3594 0 : vector<PseudoJet> pieces;
3595 0 : pieces.push_back(j1);
3596 0 : pieces.push_back(j2);
3597 0 : return join(pieces);
3598 0 : }
3599 : PseudoJet join(const PseudoJet & j1, const PseudoJet & j2, const PseudoJet & j3){
3600 0 : vector<PseudoJet> pieces;
3601 0 : pieces.push_back(j1);
3602 0 : pieces.push_back(j2);
3603 0 : pieces.push_back(j3);
3604 0 : return join(pieces);
3605 0 : }
3606 : PseudoJet join(const PseudoJet & j1, const PseudoJet & j2, const PseudoJet & j3, const PseudoJet & j4){
3607 0 : vector<PseudoJet> pieces;
3608 0 : pieces.push_back(j1);
3609 0 : pieces.push_back(j2);
3610 0 : pieces.push_back(j3);
3611 0 : pieces.push_back(j4);
3612 0 : return join(pieces);
3613 0 : }
3614 : FJCORE_END_NAMESPACE
3615 : using namespace std;
3616 : FJCORE_BEGIN_NAMESPACE // defined in fastjet/internal/base.hh
3617 : const ClusterSequence* PseudoJetStructureBase::associated_cluster_sequence() const{
3618 0 : return NULL;
3619 : }
3620 : const ClusterSequence * PseudoJetStructureBase::validated_cs() const{
3621 0 : throw Error("This PseudoJet structure is not associated with a valid ClusterSequence");
3622 0 : }
3623 : bool PseudoJetStructureBase::has_partner(const PseudoJet & /*reference */, PseudoJet & /*partner*/) const{
3624 0 : throw Error("This PseudoJet structure has no implementation for has_partner");
3625 0 : }
3626 : bool PseudoJetStructureBase::has_child(const PseudoJet & /*reference*/, PseudoJet & /*child*/) const{
3627 0 : throw Error("This PseudoJet structure has no implementation for has_child");
3628 0 : }
3629 : bool PseudoJetStructureBase::has_parents(const PseudoJet & /*reference*/, PseudoJet &/*parent1*/, PseudoJet &/*parent2*/) const{
3630 0 : throw Error("This PseudoJet structure has no implementation for has_parents");
3631 0 : }
3632 : bool PseudoJetStructureBase::object_in_jet(const PseudoJet & /*reference*/, const PseudoJet & /*jet*/) const{
3633 0 : throw Error("This PseudoJet structure has no implementation for is_inside");
3634 0 : }
3635 : vector<PseudoJet> PseudoJetStructureBase::constituents(const PseudoJet &/*reference*/) const{
3636 0 : throw Error("This PseudoJet structure has no implementation for constituents");
3637 0 : }
3638 : vector<PseudoJet> PseudoJetStructureBase::exclusive_subjets (const PseudoJet & /*reference*/, const double & /*dcut*/) const{
3639 0 : throw Error("This PseudoJet structure has no implementation for exclusive_subjets");
3640 0 : }
3641 : int PseudoJetStructureBase::n_exclusive_subjets(const PseudoJet & /*reference*/, const double & /*dcut*/) const{
3642 0 : throw Error("This PseudoJet structure has no implementation for n_exclusive_subjets");
3643 0 : }
3644 : vector<PseudoJet> PseudoJetStructureBase::exclusive_subjets_up_to (const PseudoJet & /*reference*/, int /*nsub*/) const{
3645 0 : throw Error("This PseudoJet structure has no implementation for exclusive_subjets");
3646 0 : }
3647 : double PseudoJetStructureBase::exclusive_subdmerge(const PseudoJet & /*reference*/, int /*nsub*/) const{
3648 0 : throw Error("This PseudoJet structure has no implementation for exclusive_submerge");
3649 0 : }
3650 : double PseudoJetStructureBase::exclusive_subdmerge_max(const PseudoJet & /*reference*/, int /*nsub*/) const{
3651 0 : throw Error("This PseudoJet structure has no implementation for exclusive_submerge_max");
3652 0 : }
3653 : std::vector<PseudoJet> PseudoJetStructureBase::pieces(const PseudoJet & /*reference*/) const{
3654 0 : throw Error("This PseudoJet structure has no implementation for pieces");
3655 0 : }
3656 : FJCORE_END_NAMESPACE
3657 : #include <sstream>
3658 : #include <algorithm>
3659 : using namespace std;
3660 : FJCORE_BEGIN_NAMESPACE // defined in fastjet/internal/base.hh
3661 : std::vector<PseudoJet> Selector::operator()(const std::vector<PseudoJet> & jets) const {
3662 0 : std::vector<PseudoJet> result;
3663 0 : const SelectorWorker * worker_local = validated_worker();
3664 0 : if (worker_local->applies_jet_by_jet()) {
3665 0 : for (std::vector<PseudoJet>::const_iterator jet = jets.begin();
3666 0 : jet != jets.end(); jet++) {
3667 0 : if (worker_local->pass(*jet)) result.push_back(*jet);
3668 : }
3669 0 : } else {
3670 0 : std::vector<const PseudoJet *> jetptrs(jets.size());
3671 0 : for (unsigned i = 0; i < jets.size(); i++) {
3672 0 : jetptrs[i] = & jets[i];
3673 : }
3674 0 : worker_local->terminator(jetptrs);
3675 0 : for (unsigned i = 0; i < jetptrs.size(); i++) {
3676 0 : if (jetptrs[i]) result.push_back(jets[i]);
3677 : }
3678 0 : }
3679 : return result;
3680 0 : }
3681 : unsigned int Selector::count(const std::vector<PseudoJet> & jets) const {
3682 : unsigned n = 0;
3683 0 : const SelectorWorker * worker_local = validated_worker();
3684 0 : if (worker_local->applies_jet_by_jet()) {
3685 0 : for (unsigned i = 0; i < jets.size(); i++) {
3686 0 : if (worker_local->pass(jets[i])) n++;
3687 : }
3688 0 : } else {
3689 0 : std::vector<const PseudoJet *> jetptrs(jets.size());
3690 0 : for (unsigned i = 0; i < jets.size(); i++) {
3691 0 : jetptrs[i] = & jets[i];
3692 : }
3693 0 : worker_local->terminator(jetptrs);
3694 0 : for (unsigned i = 0; i < jetptrs.size(); i++) {
3695 0 : if (jetptrs[i]) n++;
3696 : }
3697 0 : }
3698 0 : return n;
3699 0 : }
3700 : void Selector::sift(const std::vector<PseudoJet> & jets,
3701 : std::vector<PseudoJet> & jets_that_pass,
3702 : std::vector<PseudoJet> & jets_that_fail
3703 : ) const {
3704 0 : const SelectorWorker * worker_local = validated_worker();
3705 0 : jets_that_pass.clear();
3706 0 : jets_that_fail.clear();
3707 0 : if (worker_local->applies_jet_by_jet()) {
3708 0 : for (unsigned i = 0; i < jets.size(); i++) {
3709 0 : if (worker_local->pass(jets[i])) {
3710 0 : jets_that_pass.push_back(jets[i]);
3711 0 : } else {
3712 0 : jets_that_fail.push_back(jets[i]);
3713 : }
3714 : }
3715 0 : } else {
3716 0 : std::vector<const PseudoJet *> jetptrs(jets.size());
3717 0 : for (unsigned i = 0; i < jets.size(); i++) {
3718 0 : jetptrs[i] = & jets[i];
3719 : }
3720 0 : worker_local->terminator(jetptrs);
3721 0 : for (unsigned i = 0; i < jetptrs.size(); i++) {
3722 0 : if (jetptrs[i]) {
3723 0 : jets_that_pass.push_back(jets[i]);
3724 : } else {
3725 0 : jets_that_fail.push_back(jets[i]);
3726 : }
3727 : }
3728 0 : }
3729 0 : }
3730 : bool SelectorWorker::has_finite_area() const {
3731 0 : if (! is_geometric()) return false;
3732 0 : double rapmin, rapmax;
3733 0 : get_rapidity_extent(rapmin, rapmax);
3734 0 : return (rapmax != std::numeric_limits<double>::infinity())
3735 0 : && (-rapmin != std::numeric_limits<double>::infinity());
3736 0 : }
3737 0 : class SW_Identity : public SelectorWorker {
3738 : public:
3739 0 : SW_Identity(){}
3740 : virtual bool pass(const PseudoJet &) const {
3741 0 : return true;
3742 : }
3743 : virtual void terminator(vector<const PseudoJet *> &) const {
3744 0 : return;
3745 : }
3746 0 : virtual string description() const { return "Identity";}
3747 0 : virtual bool is_geometric() const { return true;}
3748 : };
3749 : Selector SelectorIdentity() {
3750 0 : return Selector(new SW_Identity);
3751 0 : }
3752 0 : class SW_Not : public SelectorWorker {
3753 : public:
3754 0 : SW_Not(const Selector & s) : _s(s) {}
3755 0 : virtual SelectorWorker* copy(){ return new SW_Not(*this);}
3756 : virtual bool pass(const PseudoJet & jet) const {
3757 0 : if (!applies_jet_by_jet())
3758 0 : throw Error("Cannot apply this selector worker to an individual jet");
3759 0 : return ! _s.pass(jet);
3760 0 : }
3761 0 : virtual bool applies_jet_by_jet() const {return _s.applies_jet_by_jet();}
3762 : virtual void terminator(vector<const PseudoJet *> & jets) const {
3763 0 : if (applies_jet_by_jet()){
3764 0 : SelectorWorker::terminator(jets);
3765 0 : return;
3766 : }
3767 0 : vector<const PseudoJet *> s_jets = jets;
3768 0 : _s.worker()->terminator(s_jets);
3769 0 : for (unsigned int i=0; i<s_jets.size(); i++){
3770 0 : if (s_jets[i]) jets[i] = NULL;
3771 : }
3772 0 : }
3773 : virtual string description() const {
3774 0 : ostringstream ostr;
3775 0 : ostr << "!(" << _s.description() << ")";
3776 0 : return ostr.str();
3777 0 : }
3778 0 : virtual bool is_geometric() const { return _s.is_geometric();}
3779 0 : virtual bool takes_reference() const { return _s.takes_reference();}
3780 0 : virtual void set_reference(const PseudoJet &ref) { _s.set_reference(ref);}
3781 : protected:
3782 : Selector _s;
3783 : };
3784 : Selector operator!(const Selector & s) {
3785 0 : return Selector(new SW_Not(s));
3786 0 : }
3787 0 : class SW_BinaryOperator: public SelectorWorker {
3788 : public:
3789 0 : SW_BinaryOperator(const Selector & s1, const Selector & s2) : _s1(s1), _s2(s2) {
3790 0 : _applies_jet_by_jet = _s1.applies_jet_by_jet() && _s2.applies_jet_by_jet();
3791 0 : _takes_reference = _s1.takes_reference() || _s2.takes_reference();
3792 0 : _is_geometric = _s1.is_geometric() && _s2.is_geometric();
3793 0 : }
3794 0 : virtual bool applies_jet_by_jet() const {return _applies_jet_by_jet;}
3795 : virtual bool takes_reference() const{
3796 0 : return _takes_reference;
3797 : }
3798 : virtual void set_reference(const PseudoJet ¢re){
3799 0 : _s1.set_reference(centre);
3800 0 : _s2.set_reference(centre);
3801 0 : }
3802 0 : virtual bool is_geometric() const { return _is_geometric;}
3803 : protected:
3804 : Selector _s1, _s2;
3805 : bool _applies_jet_by_jet;
3806 : bool _takes_reference;
3807 : bool _is_geometric;
3808 : };
3809 0 : class SW_And: public SW_BinaryOperator {
3810 : public:
3811 0 : SW_And(const Selector & s1, const Selector & s2) : SW_BinaryOperator(s1,s2){}
3812 0 : virtual SelectorWorker* copy(){ return new SW_And(*this);}
3813 : virtual bool pass(const PseudoJet & jet) const {
3814 0 : if (!applies_jet_by_jet())
3815 0 : throw Error("Cannot apply this selector worker to an individual jet");
3816 0 : return _s1.pass(jet) && _s2.pass(jet);
3817 0 : }
3818 : virtual void terminator(vector<const PseudoJet *> & jets) const {
3819 0 : if (applies_jet_by_jet()){
3820 0 : SelectorWorker::terminator(jets);
3821 0 : return;
3822 : }
3823 0 : vector<const PseudoJet *> s1_jets = jets;
3824 0 : _s1.worker()->terminator(s1_jets);
3825 0 : _s2.worker()->terminator(jets);
3826 0 : for (unsigned int i=0; i<jets.size(); i++){
3827 0 : if (! s1_jets[i]) jets[i] = NULL;
3828 : }
3829 0 : }
3830 : virtual void get_rapidity_extent(double & rapmin, double & rapmax) const {
3831 0 : double s1min, s1max, s2min, s2max;
3832 0 : _s1.get_rapidity_extent(s1min, s1max);
3833 0 : _s2.get_rapidity_extent(s2min, s2max);
3834 0 : rapmax = min(s1max, s2max);
3835 0 : rapmin = max(s1min, s2min);
3836 0 : }
3837 : virtual string description() const {
3838 0 : ostringstream ostr;
3839 0 : ostr << "(" << _s1.description() << " && " << _s2.description() << ")";
3840 0 : return ostr.str();
3841 0 : }
3842 : };
3843 : Selector operator&&(const Selector & s1, const Selector & s2) {
3844 0 : return Selector(new SW_And(s1,s2));
3845 0 : }
3846 0 : class SW_Or: public SW_BinaryOperator {
3847 : public:
3848 0 : SW_Or(const Selector & s1, const Selector & s2) : SW_BinaryOperator(s1,s2) {}
3849 0 : virtual SelectorWorker* copy(){ return new SW_Or(*this);}
3850 : virtual bool pass(const PseudoJet & jet) const {
3851 0 : if (!applies_jet_by_jet())
3852 0 : throw Error("Cannot apply this selector worker to an individual jet");
3853 0 : return _s1.pass(jet) || _s2.pass(jet);
3854 0 : }
3855 : virtual bool applies_jet_by_jet() const {
3856 0 : return _s1.applies_jet_by_jet() && _s2.applies_jet_by_jet();
3857 : }
3858 : virtual void terminator(vector<const PseudoJet *> & jets) const {
3859 0 : if (applies_jet_by_jet()){
3860 0 : SelectorWorker::terminator(jets);
3861 0 : return;
3862 : }
3863 0 : vector<const PseudoJet *> s1_jets = jets;
3864 0 : _s1.worker()->terminator(s1_jets);
3865 0 : _s2.worker()->terminator(jets);
3866 0 : for (unsigned int i=0; i<jets.size(); i++){
3867 0 : if (s1_jets[i]) jets[i] = s1_jets[i];
3868 : }
3869 0 : }
3870 : virtual string description() const {
3871 0 : ostringstream ostr;
3872 0 : ostr << "(" << _s1.description() << " || " << _s2.description() << ")";
3873 0 : return ostr.str();
3874 0 : }
3875 : virtual void get_rapidity_extent(double & rapmin, double & rapmax) const {
3876 0 : double s1min, s1max, s2min, s2max;
3877 0 : _s1.get_rapidity_extent(s1min, s1max);
3878 0 : _s2.get_rapidity_extent(s2min, s2max);
3879 0 : rapmax = max(s1max, s2max);
3880 0 : rapmin = min(s1min, s2min);
3881 0 : }
3882 : };
3883 : Selector operator ||(const Selector & s1, const Selector & s2) {
3884 0 : return Selector(new SW_Or(s1,s2));
3885 0 : }
3886 0 : class SW_Mult: public SW_And {
3887 : public:
3888 0 : SW_Mult(const Selector & s1, const Selector & s2) : SW_And(s1,s2) {}
3889 0 : virtual SelectorWorker* copy(){ return new SW_Mult(*this);}
3890 : virtual void terminator(vector<const PseudoJet *> & jets) const {
3891 0 : if (applies_jet_by_jet()){
3892 0 : SelectorWorker::terminator(jets);
3893 0 : return;
3894 : }
3895 0 : _s2.worker()->terminator(jets);
3896 0 : _s1.worker()->terminator(jets);
3897 0 : }
3898 : virtual string description() const {
3899 0 : ostringstream ostr;
3900 0 : ostr << "(" << _s1.description() << " * " << _s2.description() << ")";
3901 0 : return ostr.str();
3902 0 : }
3903 : };
3904 : Selector operator*(const Selector & s1, const Selector & s2) {
3905 0 : return Selector(new SW_Mult(s1,s2));
3906 0 : }
3907 : class QuantityBase{
3908 : public:
3909 0 : QuantityBase(double q) : _q(q){}
3910 0 : virtual ~QuantityBase(){}
3911 : virtual double operator()(const PseudoJet & jet ) const =0;
3912 : virtual string description() const =0;
3913 0 : virtual bool is_geometric() const { return false;}
3914 0 : virtual double comparison_value() const {return _q;}
3915 0 : virtual double description_value() const {return comparison_value();}
3916 : protected:
3917 : double _q;
3918 : };
3919 0 : class QuantitySquareBase : public QuantityBase{
3920 : public:
3921 0 : QuantitySquareBase(double sqrtq) : QuantityBase(sqrtq*sqrtq), _sqrtq(sqrtq){}
3922 0 : virtual double description_value() const {return _sqrtq;}
3923 : protected:
3924 : double _sqrtq;
3925 : };
3926 : template<typename QuantityType>
3927 0 : class SW_QuantityMin : public SelectorWorker{
3928 : public:
3929 0 : SW_QuantityMin(double qmin) : _qmin(qmin) {}
3930 0 : virtual bool pass(const PseudoJet & jet) const {return _qmin(jet) >= _qmin.comparison_value();}
3931 : virtual string description() const {
3932 0 : ostringstream ostr;
3933 0 : ostr << _qmin.description() << " >= " << _qmin.description_value();
3934 0 : return ostr.str();
3935 0 : }
3936 0 : virtual bool is_geometric() const { return _qmin.is_geometric();}
3937 : protected:
3938 : QuantityType _qmin; ///< the cut
3939 : };
3940 : template<typename QuantityType>
3941 0 : class SW_QuantityMax : public SelectorWorker {
3942 : public:
3943 0 : SW_QuantityMax(double qmax) : _qmax(qmax) {}
3944 0 : virtual bool pass(const PseudoJet & jet) const {return _qmax(jet) <= _qmax.comparison_value();}
3945 : virtual string description() const {
3946 0 : ostringstream ostr;
3947 0 : ostr << _qmax.description() << " <= " << _qmax.description_value();
3948 0 : return ostr.str();
3949 0 : }
3950 0 : virtual bool is_geometric() const { return _qmax.is_geometric();}
3951 : protected:
3952 : QuantityType _qmax; ///< the cut
3953 : };
3954 : template<typename QuantityType>
3955 0 : class SW_QuantityRange : public SelectorWorker {
3956 : public:
3957 0 : SW_QuantityRange(double qmin, double qmax) : _qmin(qmin), _qmax(qmax) {}
3958 : virtual bool pass(const PseudoJet & jet) const {
3959 0 : double q = _qmin(jet); // we could identically use _qmax
3960 0 : return (q >= _qmin.comparison_value()) && (q <= _qmax.comparison_value());
3961 : }
3962 : virtual string description() const {
3963 0 : ostringstream ostr;
3964 0 : ostr << _qmin.description_value() << " <= " << _qmin.description() << " <= " << _qmax.description_value();
3965 0 : return ostr.str();
3966 0 : }
3967 0 : virtual bool is_geometric() const { return _qmin.is_geometric();}
3968 : protected:
3969 : QuantityType _qmin; // the lower cut
3970 : QuantityType _qmax; // the upper cut
3971 : };
3972 0 : class QuantityPt2 : public QuantitySquareBase{
3973 : public:
3974 0 : QuantityPt2(double pt) : QuantitySquareBase(pt){}
3975 0 : virtual double operator()(const PseudoJet & jet ) const { return jet.perp2();}
3976 0 : virtual string description() const {return "pt";}
3977 : };
3978 : Selector SelectorPtMin(double ptmin) {
3979 0 : return Selector(new SW_QuantityMin<QuantityPt2>(ptmin));
3980 0 : }
3981 : Selector SelectorPtMax(double ptmax) {
3982 0 : return Selector(new SW_QuantityMax<QuantityPt2>(ptmax));
3983 0 : }
3984 : Selector SelectorPtRange(double ptmin, double ptmax) {
3985 0 : return Selector(new SW_QuantityRange<QuantityPt2>(ptmin, ptmax));
3986 0 : }
3987 0 : class QuantityEt2 : public QuantitySquareBase{
3988 : public:
3989 0 : QuantityEt2(double Et) : QuantitySquareBase(Et){}
3990 0 : virtual double operator()(const PseudoJet & jet ) const { return jet.Et2();}
3991 0 : virtual string description() const {return "Et";}
3992 : };
3993 : Selector SelectorEtMin(double Etmin) {
3994 0 : return Selector(new SW_QuantityMin<QuantityEt2>(Etmin));
3995 0 : }
3996 : Selector SelectorEtMax(double Etmax) {
3997 0 : return Selector(new SW_QuantityMax<QuantityEt2>(Etmax));
3998 0 : }
3999 : Selector SelectorEtRange(double Etmin, double Etmax) {
4000 0 : return Selector(new SW_QuantityRange<QuantityEt2>(Etmin, Etmax));
4001 0 : }
4002 0 : class QuantityE : public QuantityBase{
4003 : public:
4004 0 : QuantityE(double E) : QuantityBase(E){}
4005 0 : virtual double operator()(const PseudoJet & jet ) const { return jet.E();}
4006 0 : virtual string description() const {return "E";}
4007 : };
4008 : Selector SelectorEMin(double Emin) {
4009 0 : return Selector(new SW_QuantityMin<QuantityE>(Emin));
4010 0 : }
4011 : Selector SelectorEMax(double Emax) {
4012 0 : return Selector(new SW_QuantityMax<QuantityE>(Emax));
4013 0 : }
4014 : Selector SelectorERange(double Emin, double Emax) {
4015 0 : return Selector(new SW_QuantityRange<QuantityE>(Emin, Emax));
4016 0 : }
4017 0 : class QuantityM2 : public QuantitySquareBase{
4018 : public:
4019 0 : QuantityM2(double m) : QuantitySquareBase(m){}
4020 0 : virtual double operator()(const PseudoJet & jet ) const { return jet.m2();}
4021 0 : virtual string description() const {return "mass";}
4022 : };
4023 : Selector SelectorMassMin(double mmin) {
4024 0 : return Selector(new SW_QuantityMin<QuantityM2>(mmin));
4025 0 : }
4026 : Selector SelectorMassMax(double mmax) {
4027 0 : return Selector(new SW_QuantityMax<QuantityM2>(mmax));
4028 0 : }
4029 : Selector SelectorMassRange(double mmin, double mmax) {
4030 0 : return Selector(new SW_QuantityRange<QuantityM2>(mmin, mmax));
4031 0 : }
4032 0 : class QuantityRap : public QuantityBase{
4033 : public:
4034 0 : QuantityRap(double rap) : QuantityBase(rap){}
4035 0 : virtual double operator()(const PseudoJet & jet ) const { return jet.rap();}
4036 0 : virtual string description() const {return "rap";}
4037 0 : virtual bool is_geometric() const { return true;}
4038 : };
4039 0 : class SW_RapMin : public SW_QuantityMin<QuantityRap>{
4040 : public:
4041 0 : SW_RapMin(double rapmin) : SW_QuantityMin<QuantityRap>(rapmin){}
4042 : virtual void get_rapidity_extent(double &rapmin, double & rapmax) const{
4043 0 : rapmax = std::numeric_limits<double>::max();
4044 0 : rapmin = _qmin.comparison_value();
4045 0 : }
4046 : };
4047 0 : class SW_RapMax : public SW_QuantityMax<QuantityRap>{
4048 : public:
4049 0 : SW_RapMax(double rapmax) : SW_QuantityMax<QuantityRap>(rapmax){}
4050 : virtual void get_rapidity_extent(double &rapmin, double & rapmax) const{
4051 0 : rapmax = _qmax.comparison_value();
4052 0 : rapmin = -std::numeric_limits<double>::max();
4053 0 : }
4054 : };
4055 0 : class SW_RapRange : public SW_QuantityRange<QuantityRap>{
4056 : public:
4057 0 : SW_RapRange(double rapmin, double rapmax) : SW_QuantityRange<QuantityRap>(rapmin, rapmax){
4058 0 : assert(rapmin<=rapmax);
4059 0 : }
4060 : virtual void get_rapidity_extent(double &rapmin, double & rapmax) const{
4061 0 : rapmax = _qmax.comparison_value();
4062 0 : rapmin = _qmin.comparison_value();
4063 0 : }
4064 0 : virtual bool has_known_area() const { return true;} ///< the area is analytically known
4065 : virtual double known_area() const {
4066 0 : return twopi * (_qmax.comparison_value()-_qmin.comparison_value());
4067 : }
4068 : };
4069 : Selector SelectorRapMin(double rapmin) {
4070 0 : return Selector(new SW_RapMin(rapmin));
4071 0 : }
4072 : Selector SelectorRapMax(double rapmax) {
4073 0 : return Selector(new SW_RapMax(rapmax));
4074 0 : }
4075 : Selector SelectorRapRange(double rapmin, double rapmax) {
4076 0 : return Selector(new SW_RapRange(rapmin, rapmax));
4077 0 : }
4078 0 : class QuantityAbsRap : public QuantityBase{
4079 : public:
4080 0 : QuantityAbsRap(double absrap) : QuantityBase(absrap){}
4081 0 : virtual double operator()(const PseudoJet & jet ) const { return abs(jet.rap());}
4082 0 : virtual string description() const {return "|rap|";}
4083 0 : virtual bool is_geometric() const { return true;}
4084 : };
4085 0 : class SW_AbsRapMax : public SW_QuantityMax<QuantityAbsRap>{
4086 : public:
4087 0 : SW_AbsRapMax(double absrapmax) : SW_QuantityMax<QuantityAbsRap>(absrapmax){}
4088 : virtual void get_rapidity_extent(double &rapmin, double & rapmax) const{
4089 0 : rapmax = _qmax.comparison_value();
4090 0 : rapmin = -_qmax.comparison_value();
4091 0 : }
4092 0 : virtual bool has_known_area() const { return true;} ///< the area is analytically known
4093 : virtual double known_area() const {
4094 0 : return twopi * 2 * _qmax.comparison_value();
4095 : }
4096 : };
4097 0 : class SW_AbsRapRange : public SW_QuantityRange<QuantityAbsRap>{
4098 : public:
4099 0 : SW_AbsRapRange(double absrapmin, double absrapmax) : SW_QuantityRange<QuantityAbsRap>(absrapmin, absrapmax){}
4100 : virtual void get_rapidity_extent(double &rapmin, double & rapmax) const{
4101 0 : rapmax = _qmax.comparison_value();
4102 0 : rapmin = -_qmax.comparison_value();
4103 0 : }
4104 0 : virtual bool has_known_area() const { return true;} ///< the area is analytically known
4105 : virtual double known_area() const {
4106 0 : return twopi * 2 * (_qmax.comparison_value()-max(_qmin.comparison_value(),0.0)); // this should handle properly absrapmin<0
4107 : }
4108 : };
4109 : Selector SelectorAbsRapMin(double absrapmin) {
4110 0 : return Selector(new SW_QuantityMin<QuantityAbsRap>(absrapmin));
4111 0 : }
4112 : Selector SelectorAbsRapMax(double absrapmax) {
4113 0 : return Selector(new SW_AbsRapMax(absrapmax));
4114 0 : }
4115 : Selector SelectorAbsRapRange(double rapmin, double rapmax) {
4116 0 : return Selector(new SW_AbsRapRange(rapmin, rapmax));
4117 0 : }
4118 0 : class QuantityEta : public QuantityBase{
4119 : public:
4120 0 : QuantityEta(double eta) : QuantityBase(eta){}
4121 0 : virtual double operator()(const PseudoJet & jet ) const { return jet.eta();}
4122 0 : virtual string description() const {return "eta";}
4123 : };
4124 : Selector SelectorEtaMin(double etamin) {
4125 0 : return Selector(new SW_QuantityMin<QuantityEta>(etamin));
4126 0 : }
4127 : Selector SelectorEtaMax(double etamax) {
4128 0 : return Selector(new SW_QuantityMax<QuantityEta>(etamax));
4129 0 : }
4130 : Selector SelectorEtaRange(double etamin, double etamax) {
4131 0 : return Selector(new SW_QuantityRange<QuantityEta>(etamin, etamax));
4132 0 : }
4133 0 : class QuantityAbsEta : public QuantityBase{
4134 : public:
4135 0 : QuantityAbsEta(double abseta) : QuantityBase(abseta){}
4136 0 : virtual double operator()(const PseudoJet & jet ) const { return abs(jet.eta());}
4137 0 : virtual string description() const {return "|eta|";}
4138 0 : virtual bool is_geometric() const { return true;}
4139 : };
4140 : Selector SelectorAbsEtaMin(double absetamin) {
4141 0 : return Selector(new SW_QuantityMin<QuantityAbsEta>(absetamin));
4142 0 : }
4143 : Selector SelectorAbsEtaMax(double absetamax) {
4144 0 : return Selector(new SW_QuantityMax<QuantityAbsEta>(absetamax));
4145 0 : }
4146 : Selector SelectorAbsEtaRange(double absetamin, double absetamax) {
4147 0 : return Selector(new SW_QuantityRange<QuantityAbsEta>(absetamin, absetamax));
4148 0 : }
4149 0 : class SW_PhiRange : public SelectorWorker {
4150 : public:
4151 0 : SW_PhiRange(double phimin, double phimax) : _phimin(phimin), _phimax(phimax){
4152 0 : assert(_phimin<_phimax);
4153 0 : assert(_phimin>-twopi);
4154 0 : assert(_phimax<2*twopi);
4155 0 : _phispan = _phimax - _phimin;
4156 0 : }
4157 : virtual bool pass(const PseudoJet & jet) const {
4158 0 : double dphi=jet.phi()-_phimin;
4159 0 : if (dphi >= twopi) dphi -= twopi;
4160 0 : if (dphi < 0) dphi += twopi;
4161 0 : return (dphi <= _phispan);
4162 : }
4163 : virtual string description() const {
4164 0 : ostringstream ostr;
4165 0 : ostr << _phimin << " <= phi <= " << _phimax;
4166 0 : return ostr.str();
4167 0 : }
4168 0 : virtual bool is_geometric() const { return true;}
4169 : protected:
4170 : double _phimin; // the lower cut
4171 : double _phimax; // the upper cut
4172 : double _phispan; // the span of the range
4173 : };
4174 : Selector SelectorPhiRange(double phimin, double phimax) {
4175 0 : return Selector(new SW_PhiRange(phimin, phimax));
4176 0 : }
4177 0 : class SW_RapPhiRange : public SW_And{
4178 : public:
4179 : SW_RapPhiRange(double rapmin, double rapmax, double phimin, double phimax)
4180 0 : : SW_And(SelectorRapRange(rapmin, rapmax), SelectorPhiRange(phimin, phimax)){
4181 0 : _known_area = ((phimax-phimin > twopi) ? twopi : phimax-phimin) * (rapmax-rapmin);
4182 0 : }
4183 : virtual double known_area() const{
4184 0 : return _known_area;
4185 : }
4186 : protected:
4187 : double _known_area;
4188 : };
4189 : Selector SelectorRapPhiRange(double rapmin, double rapmax, double phimin, double phimax) {
4190 0 : return Selector(new SW_RapPhiRange(rapmin, rapmax, phimin, phimax));
4191 0 : }
4192 0 : class SW_NHardest : public SelectorWorker {
4193 : public:
4194 0 : SW_NHardest(unsigned int n) : _n(n) {};
4195 : virtual bool pass(const PseudoJet &) const {
4196 0 : if (!applies_jet_by_jet())
4197 0 : throw Error("Cannot apply this selector worker to an individual jet");
4198 0 : return false;
4199 0 : }
4200 : virtual void terminator(vector<const PseudoJet *> & jets) const {
4201 0 : if (jets.size() < _n) return;
4202 0 : vector<double> minus_pt2(jets.size());
4203 0 : vector<unsigned int> indices(jets.size());
4204 0 : for (unsigned int i=0; i<jets.size(); i++){
4205 0 : indices[i] = i;
4206 0 : minus_pt2[i] = jets[i] ? -jets[i]->perp2() : 0.0;
4207 : }
4208 0 : IndexedSortHelper sort_helper(& minus_pt2);
4209 0 : partial_sort(indices.begin(), indices.begin()+_n, indices.end(), sort_helper);
4210 0 : for (unsigned int i=_n; i<jets.size(); i++)
4211 0 : jets[indices[i]] = NULL;
4212 0 : }
4213 0 : virtual bool applies_jet_by_jet() const {return false;}
4214 : virtual string description() const {
4215 0 : ostringstream ostr;
4216 0 : ostr << _n << " hardest";
4217 0 : return ostr.str();
4218 0 : }
4219 : protected:
4220 : unsigned int _n;
4221 : };
4222 : Selector SelectorNHardest(unsigned int n) {
4223 0 : return Selector(new SW_NHardest(n));
4224 0 : }
4225 0 : class SW_WithReference : public SelectorWorker{
4226 : public:
4227 0 : SW_WithReference() : _is_initialised(false){};
4228 0 : virtual bool takes_reference() const { return true;}
4229 : virtual void set_reference(const PseudoJet ¢re){
4230 0 : _is_initialised = true;
4231 0 : _reference = centre;
4232 0 : }
4233 : protected:
4234 : PseudoJet _reference;
4235 : bool _is_initialised;
4236 : };
4237 0 : class SW_Circle : public SW_WithReference {
4238 : public:
4239 0 : SW_Circle(const double &radius) : _radius2(radius*radius) {}
4240 0 : virtual SelectorWorker* copy(){ return new SW_Circle(*this);}
4241 : virtual bool pass(const PseudoJet & jet) const {
4242 0 : if (! _is_initialised)
4243 0 : throw Error("To use a SelectorCircle (or any selector that requires a reference), you first have to call set_reference(...)");
4244 0 : return jet.squared_distance(_reference) <= _radius2;
4245 0 : }
4246 : virtual string description() const {
4247 0 : ostringstream ostr;
4248 0 : ostr << "distance from the centre <= " << sqrt(_radius2);
4249 0 : return ostr.str();
4250 0 : }
4251 : virtual void get_rapidity_extent(double & rapmin, double & rapmax) const{
4252 0 : if (! _is_initialised)
4253 0 : throw Error("To use a SelectorCircle (or any selector that requires a reference), you first have to call set_reference(...)");
4254 0 : rapmax = _reference.rap()+sqrt(_radius2);
4255 0 : rapmin = _reference.rap()-sqrt(_radius2);
4256 0 : }
4257 0 : virtual bool is_geometric() const { return true;} ///< implies a finite area
4258 0 : virtual bool has_finite_area() const { return true;} ///< regardless of the reference
4259 0 : virtual bool has_known_area() const { return true;} ///< the area is analytically known
4260 : virtual double known_area() const {
4261 0 : return pi * _radius2;
4262 : }
4263 : protected:
4264 : double _radius2;
4265 : };
4266 : Selector SelectorCircle(const double & radius) {
4267 0 : return Selector(new SW_Circle(radius));
4268 0 : }
4269 0 : class SW_Doughnut : public SW_WithReference {
4270 : public:
4271 0 : SW_Doughnut(const double &radius_in, const double &radius_out)
4272 0 : : _radius_in2(radius_in*radius_in), _radius_out2(radius_out*radius_out) {}
4273 0 : virtual SelectorWorker* copy(){ return new SW_Doughnut(*this);}
4274 : virtual bool pass(const PseudoJet & jet) const {
4275 0 : if (! _is_initialised)
4276 0 : throw Error("To use a SelectorDoughnut (or any selector that requires a reference), you first have to call set_reference(...)");
4277 0 : double distance2 = jet.squared_distance(_reference);
4278 0 : return (distance2 <= _radius_out2) && (distance2 >= _radius_in2);
4279 0 : }
4280 : virtual string description() const {
4281 0 : ostringstream ostr;
4282 0 : ostr << sqrt(_radius_in2) << " <= distance from the centre <= " << sqrt(_radius_out2);
4283 0 : return ostr.str();
4284 0 : }
4285 : virtual void get_rapidity_extent(double & rapmin, double & rapmax) const{
4286 0 : if (! _is_initialised)
4287 0 : throw Error("To use a SelectorDoughnut (or any selector that requires a reference), you first have to call set_reference(...)");
4288 0 : rapmax = _reference.rap()+sqrt(_radius_out2);
4289 0 : rapmin = _reference.rap()-sqrt(_radius_out2);
4290 0 : }
4291 0 : virtual bool is_geometric() const { return true;} ///< implies a finite area
4292 0 : virtual bool has_finite_area() const { return true;} ///< regardless of the reference
4293 0 : virtual bool has_known_area() const { return true;} ///< the area is analytically known
4294 : virtual double known_area() const {
4295 0 : return pi * (_radius_out2-_radius_in2);
4296 : }
4297 : protected:
4298 : double _radius_in2, _radius_out2;
4299 : };
4300 : Selector SelectorDoughnut(const double & radius_in, const double & radius_out) {
4301 0 : return Selector(new SW_Doughnut(radius_in, radius_out));
4302 0 : }
4303 0 : class SW_Strip : public SW_WithReference {
4304 : public:
4305 0 : SW_Strip(const double &delta) : _delta(delta) {}
4306 0 : virtual SelectorWorker* copy(){ return new SW_Strip(*this);}
4307 : virtual bool pass(const PseudoJet & jet) const {
4308 0 : if (! _is_initialised)
4309 0 : throw Error("To use a SelectorStrip (or any selector that requires a reference), you first have to call set_reference(...)");
4310 0 : return abs(jet.rap()-_reference.rap()) <= _delta;
4311 0 : }
4312 : virtual string description() const {
4313 0 : ostringstream ostr;
4314 0 : ostr << "|rap - rap_reference| <= " << _delta;
4315 0 : return ostr.str();
4316 0 : }
4317 : virtual void get_rapidity_extent(double & rapmin, double & rapmax) const{
4318 0 : if (! _is_initialised)
4319 0 : throw Error("To use a SelectorStrip (or any selector that requires a reference), you first have to call set_reference(...)");
4320 0 : rapmax = _reference.rap()+_delta;
4321 0 : rapmin = _reference.rap()-_delta;
4322 0 : }
4323 0 : virtual bool is_geometric() const { return true;} ///< implies a finite area
4324 0 : virtual bool has_finite_area() const { return true;} ///< regardless of the reference
4325 0 : virtual bool has_known_area() const { return true;} ///< the area is analytically known
4326 : virtual double known_area() const {
4327 0 : return twopi * 2 * _delta;
4328 : }
4329 : protected:
4330 : double _delta;
4331 : };
4332 : Selector SelectorStrip(const double & half_width) {
4333 0 : return Selector(new SW_Strip(half_width));
4334 0 : }
4335 0 : class SW_Rectangle : public SW_WithReference {
4336 : public:
4337 0 : SW_Rectangle(const double &delta_rap, const double &delta_phi)
4338 0 : : _delta_rap(delta_rap), _delta_phi(delta_phi) {}
4339 0 : virtual SelectorWorker* copy(){ return new SW_Rectangle(*this);}
4340 : virtual bool pass(const PseudoJet & jet) const {
4341 0 : if (! _is_initialised)
4342 0 : throw Error("To use a SelectorRectangle (or any selector that requires a reference), you first have to call set_reference(...)");
4343 0 : return (abs(jet.rap()-_reference.rap()) <= _delta_rap) && (abs(jet.delta_phi_to(_reference)) <= _delta_phi);
4344 0 : }
4345 : virtual string description() const {
4346 0 : ostringstream ostr;
4347 0 : ostr << "|rap - rap_reference| <= " << _delta_rap << " && |phi - phi_reference| <= " << _delta_phi ;
4348 0 : return ostr.str();
4349 0 : }
4350 : virtual void get_rapidity_extent(double & rapmin, double & rapmax) const{
4351 0 : if (! _is_initialised)
4352 0 : throw Error("To use a SelectorRectangle (or any selector that requires a reference), you first have to call set_reference(...)");
4353 0 : rapmax = _reference.rap()+_delta_rap;
4354 0 : rapmin = _reference.rap()-_delta_rap;
4355 0 : }
4356 0 : virtual bool is_geometric() const { return true;} ///< implies a finite area
4357 0 : virtual bool has_finite_area() const { return true;} ///< regardless of the reference
4358 0 : virtual bool has_known_area() const { return true;} ///< the area is analytically known
4359 : virtual double known_area() const {
4360 0 : return 4 * _delta_rap * _delta_phi;
4361 : }
4362 : protected:
4363 : double _delta_rap, _delta_phi;
4364 : };
4365 : Selector SelectorRectangle(const double & half_rap_width, const double & half_phi_width) {
4366 0 : return Selector(new SW_Rectangle(half_rap_width, half_phi_width));
4367 0 : }
4368 0 : class SW_PtFractionMin : public SW_WithReference {
4369 : public:
4370 0 : SW_PtFractionMin(double fraction) : _fraction2(fraction*fraction){}
4371 0 : virtual SelectorWorker* copy(){ return new SW_PtFractionMin(*this);}
4372 : virtual bool pass(const PseudoJet & jet) const {
4373 0 : if (! _is_initialised)
4374 0 : throw Error("To use a SelectorPtFractionMin (or any selector that requires a reference), you first have to call set_reference(...)");
4375 0 : return (jet.perp2() >= _fraction2*_reference.perp2());
4376 0 : }
4377 : virtual string description() const {
4378 0 : ostringstream ostr;
4379 0 : ostr << "pt >= " << sqrt(_fraction2) << "* pt_ref";
4380 0 : return ostr.str();
4381 0 : }
4382 : protected:
4383 : double _fraction2;
4384 : };
4385 : Selector SelectorPtFractionMin(double fraction){
4386 0 : return Selector(new SW_PtFractionMin(fraction));
4387 0 : }
4388 0 : class SW_IsZero : public SelectorWorker {
4389 : public:
4390 0 : SW_IsZero(){}
4391 : virtual bool pass(const PseudoJet & jet) const {
4392 0 : return jet==0;
4393 : }
4394 0 : virtual string description() const { return "zero";}
4395 : };
4396 : Selector SelectorIsZero(){
4397 0 : return Selector(new SW_IsZero());
4398 0 : }
4399 : Selector & Selector::operator &=(const Selector & b){
4400 0 : _worker.reset(new SW_And(*this, b));
4401 0 : return *this;
4402 0 : }
4403 : Selector & Selector::operator |=(const Selector & b){
4404 0 : _worker.reset(new SW_Or(*this, b));
4405 0 : return *this;
4406 0 : }
4407 : FJCORE_END_NAMESPACE // defined in fastjet/internal/base.hh
|
