79 #ifndef __FJCORE_VERSION_HH__
80 #define __FJCORE_VERSION_HH__
85 #endif // __FJCORE_VERSION_HH__
86 #ifndef __FJCORE_CLUSTERQUENCE_N2_ICC__
87 #define __FJCORE_CLUSTERQUENCE_N2_ICC__
91 BJ * briefjets =
new BJ[
n];
92 BJ * jetA = briefjets, * jetB;
93 for (
int i = 0;
i<
n;
i++) {
94 _bj_set_jetinfo(jetA,
i);
98 BJ * head = briefjets;
99 for (jetA = head + 1; jetA != tail; jetA++) {
100 _bj_set_NN_crosscheck(jetA, head, jetA);
102 double * diJ =
new double[
n];
104 for (
int i = 0;
i <
n;
i++) {
105 diJ[
i] = _bj_diJ(jetA);
108 int history_location = n-1;
109 while (tail != head) {
110 double diJ_min = diJ[0];
112 for (
int i = 1;
i <
n;
i++) {
113 if (diJ[
i] < diJ_min) {diJ_min_jet =
i; diJ_min = diJ[
i];}
116 jetA = & briefjets[diJ_min_jet];
117 jetB =
static_cast<BJ *
>(jetA->NN);
122 _do_ij_recombination_step(jetA->_jets_index, jetB->_jets_index, diJ_min, nn);
123 _bj_set_jetinfo(jetB, nn);
125 _do_iB_recombination_step(jetA->_jets_index, diJ_min);
129 diJ[jetA - head] = diJ[tail-head];
130 for (BJ * jetI = head; jetI != tail; jetI++) {
131 if (jetI->NN == jetA || jetI->NN == jetB) {
132 _bj_set_NN_nocross(jetI, head, tail);
133 diJ[jetI-head] = _bj_diJ(jetI);
136 double dist = _bj_dist(jetI,jetB);
137 if (dist < jetI->NN_dist) {
139 jetI->NN_dist =
dist;
141 diJ[jetI-head] = _bj_diJ(jetI);
144 if (dist < jetB->NN_dist) {
146 jetB->NN_dist =
dist;
150 if (jetI->NN == tail) {jetI->NN = jetA;}
152 if (jetB != NULL) {diJ[jetB-head] = _bj_diJ(jetB);}
158 #endif // __FJCORE_CLUSTERQUENCE_N2_ICC__
159 #ifndef __FJCORE_DYNAMICNEARESTNEIGHBOURS_HH__
160 #define __FJCORE_DYNAMICNEARESTNEIGHBOURS_HH__
185 virtual bool Valid(
const int index)
const = 0;
187 const std::vector<EtaPhi> & points_to_add,
188 std::vector<int> & indices_added,
189 std::vector<int> & indices_of_updated_neighbours) = 0;
191 std::vector<int> & indices_of_updated_neighbours) {
192 std::vector<int> indices_added;
193 std::vector<EtaPhi> points_to_add;
194 std::vector<int> indices_to_remove(1);
195 indices_to_remove[0] =
index;
197 indices_of_updated_neighbours
203 std::vector<int> & indices_of_updated_neighbours) {
204 std::vector<int> indices_added(1);
205 std::vector<EtaPhi> points_to_add(1);
206 std::vector<int> indices_to_remove(2);
207 indices_to_remove[0] =
index1;
208 indices_to_remove[1] =
index2;
209 points_to_add[0] = newpoint;
211 indices_of_updated_neighbours
213 index3 = indices_added[0];
218 #endif // __FJCORE_DYNAMICNEARESTNEIGHBOURS_HH__
219 #ifndef __FJCORE_SEARCHTREE_HH__
220 #define __FJCORE_SEARCHTREE_HH__
231 SearchTree(
const std::vector<T> & init,
unsigned int max_size);
232 void remove(
unsigned node_index);
242 #ifdef __FJCORE_SEARCHTREE_TRACK_DEPTH
243 inline unsigned int max_depth()
const {
return _max_depth;};
247 int loc(
const Node *
node)
const ;
260 unsigned int left_edge,
unsigned int right_edge,
262 #ifdef __FJCORE_SEARCHTREE_TRACK_DEPTH
263 unsigned int _max_depth;
270 return ((parent==0) && (
left==0) && (right==0));};
276 void reset_parents_link_to_me(
Node * XX);
285 if (
parent == NULL) {
return;}
300 _node = _node->successor;
304 _node = _node->successor;
307 _node = _node->predecessor;
311 _node = _node->predecessor;
330 _node = _node->successor;
334 _node = _node->successor;
337 _node = _node->predecessor;
341 _node = _node->predecessor;
353 unsigned int max_size) :
357 for (
unsigned int i = init.size();
i < max_size;
i++) {
363 _nodes(init.
size()), _available_nodes(0) {
369 unsigned n = init.size();
371 #ifdef __FJCORE_SEARCHTREE_TRACK_DEPTH
374 for (
unsigned int i = 1;
i<
n;
i++) {
377 for(
unsigned int i = 0;
i <
n;
i++) {
378 _nodes[
i].value = init[
i];
379 _nodes[
i].predecessor = (& (_nodes[
i])) - 1;
380 _nodes[
i].successor = (& (_nodes[
i])) + 1;
381 _nodes[
i].nullify_treelinks();
383 _nodes[0].predecessor = (& (_nodes[n-1]));
384 _nodes[n-1].successor = (& (_nodes[0]));
385 unsigned int scale = (n+1)/2;
386 unsigned int top =
std::min(n-1,scale);
387 _nodes[top].parent = NULL;
388 _top_node = &(_nodes[top]);
389 _do_initial_connections(top, scale, 0, n, 0);
392 -999 : node - &(_nodes[0]);}
394 unsigned int this_one,
396 unsigned int left_edge,
397 unsigned int right_edge,
400 #ifdef __FJCORE_SEARCHTREE_TRACK_DEPTH
401 _max_depth = max(depth, _max_depth);
403 unsigned int ref_new_scale = (scale+1)/2;
404 unsigned new_scale = ref_new_scale;
405 bool did_child =
false;
407 int left = this_one - new_scale;
408 if (left >= static_cast<int>(left_edge)
409 && _nodes[
left].treelinks_null() ) {
410 _nodes[
left].parent = &(_nodes[this_one]);
411 _nodes[this_one].left = &(_nodes[
left]);
412 _do_initial_connections(left, new_scale, left_edge, this_one, depth+1);
416 unsigned int old_new_scale = new_scale;
417 new_scale = (old_new_scale + 1)/2;
418 if (new_scale == old_new_scale)
break;
420 if (!did_child) {_nodes[this_one].left = NULL;}
421 new_scale = ref_new_scale;
424 unsigned int right = this_one + new_scale;
425 if (right < right_edge && _nodes[right].treelinks_null()) {
426 _nodes[right].parent = &(_nodes[this_one]);
427 _nodes[this_one].right = &(_nodes[right]);
428 _do_initial_connections(right, new_scale, this_one+1,right_edge,depth+1);
432 unsigned int old_new_scale = new_scale;
433 new_scale = (old_new_scale + 1)/2;
434 if (new_scale == old_new_scale)
break;
436 if (!did_child) {_nodes[this_one].right = NULL;}
439 remove(&(_nodes[node_index]));
449 if (node->
left == NULL && node->
right == NULL) {
451 }
else if (node->
left != NULL && node->
right == NULL){
454 if (_top_node == node) {_top_node = node->
left;}
455 }
else if (node->
left == NULL && node->
right != NULL){
458 if (_top_node == node) {_top_node = node->
right;}
461 bool use_predecessor = (_n_removes % 2 == 1);
462 if (use_predecessor) {
464 assert(replacement->right == NULL);
465 if (replacement != node->
left) {
466 if (replacement->left != NULL) {
467 replacement->left->parent = replacement->parent;}
468 replacement->reset_parents_link_to_me(replacement->left);
469 replacement->left = node->
left;
475 assert(replacement->left == NULL);
476 if (replacement != node->
right) {
477 if (replacement->right != NULL) {
478 replacement->right->parent = replacement->parent;}
479 replacement->reset_parents_link_to_me(replacement->right);
480 replacement->right = node->
right;
486 if (node->
left != replacement) {node->
left->
parent = replacement;}
488 if (_top_node == node) {_top_node = replacement;}
494 _available_nodes.push_back(node);
497 assert(_available_nodes.size() > 0);
498 Node * node = _available_nodes.back();
499 _available_nodes.pop_back();
501 Node * location = _top_node;
502 Node * old_location = NULL;
504 #ifdef __FJCORE_SEARCHTREE_TRACK_DEPTH
505 unsigned int depth = 0;
507 while(location != NULL) {
508 #ifdef __FJCORE_SEARCHTREE_TRACK_DEPTH
511 old_location = location;
512 on_left = value < location->
value;
513 if (on_left) {location = location->
left;}
514 else {location = location->
right;}
516 #ifdef __FJCORE_SEARCHTREE_TRACK_DEPTH
517 _max_depth = max(depth, _max_depth);
519 node->parent = old_location;
525 if (node->predecessor != NULL) {
526 node->successor = node->predecessor->successor;
527 node->predecessor->successor = node;
528 node->successor->predecessor = node;
530 node->successor = _find_successor(node);
531 assert(node->successor != NULL);
532 node->predecessor = node->successor->predecessor;
533 node->successor->predecessor = node;
534 node->predecessor->successor = node;
539 verify_structure_linear();
540 const Node * left_limit = _top_node;
541 while (left_limit->
left != NULL) {left_limit = left_limit->
left;}
542 const Node * right_limit = _top_node;
543 while (right_limit->
right != NULL) {right_limit = right_limit->
right;}
544 verify_structure_recursive(_top_node, left_limit, right_limit);
555 if (left != left_limit) {
556 verify_structure_recursive(left, left_limit, element);}
561 if (right != right_limit) {
562 verify_structure_recursive(right, element, right_limit);}
568 for(
unsigned i = 0;
i < _nodes.size();
i++) {
571 if (node->
parent == NULL) {
576 if (node->
left != NULL) {
578 if (node->
right != NULL) {
581 assert(n_top == 1 || (n_top == 0 &&
size() <= 1) );
582 assert(n_null == _available_nodes.size() ||
583 (n_null == _available_nodes.size() + 1 &&
size() == 1));
587 if (node->
left != NULL) {
588 newnode = node->
left;
589 while(newnode->
right != NULL) {newnode = newnode->
right;}
594 while(newnode != NULL) {
595 if (newnode->
right == lastnode) {
return newnode;}
597 newnode = newnode->
parent;
604 if (node->
right != NULL) {
605 newnode = node->
right;
606 while(newnode->
left != NULL) {newnode = newnode->
left;}
611 while(newnode != NULL) {
612 if (newnode->
left == lastnode) {
return newnode;}
614 newnode = newnode->
parent;
622 int n = _nodes.size();
623 for(; node - base_node <
n ; node++) {
634 #endif // __FJCORE_SEARCHTREE_HH__
635 #ifndef __FJCORE_MINHEAP__HH__
636 #define __FJCORE_MINHEAP__HH__
654 void remove(
unsigned int loc) {
655 update(
loc,std::numeric_limits<double>::max());};
656 void update(
unsigned int,
double);
665 #endif // __FJCORE_MINHEAP__HH__
666 #ifndef __FJCORE_CLOSESTPAIR2DBASE__HH__
667 #define __FJCORE_CLOSESTPAIR2DBASE__HH__
681 return Coord2D(factor*coord.
x,factor*coord.
y);
684 return Coord2D(
x / divisor,
y / divisor);};
686 double dx = a.
x - b.
x,
dy = a.
y-b.
y;
690 double dx =
x - b.
x,
dy =
y-b.
y;
696 virtual void closest_pair(
unsigned int & ID1,
unsigned int & ID2,
697 double & distance2)
const = 0;
698 virtual void remove(
unsigned int ID) = 0;
700 virtual unsigned int replace(
unsigned int ID1,
unsigned int ID2,
704 unsigned new_ID =
insert(position);
707 virtual void replace_many(
const std::vector<unsigned int> & IDs_to_remove,
708 const std::vector<Coord2D> & new_positions,
709 std::vector<unsigned int> & new_IDs) {
710 for(
unsigned i = 0;
i < IDs_to_remove.size();
i++) {
711 remove(IDs_to_remove[
i]);}
713 for(
unsigned i = 0;
i < new_positions.size();
i++) {
714 new_IDs.push_back(
insert(new_positions[
i]));}
716 virtual unsigned int size() = 0;
720 #endif // __FJCORE_CLOSESTPAIR2DBASE__HH__
721 #ifndef __FJCORE_CLOSESTPAIR2D__HH__
722 #define __FJCORE_CLOSESTPAIR2D__HH__
731 _initialize(positions, left_corner, right_corner, positions.size());
735 const unsigned int max_size) {
736 _initialize(positions, left_corner, right_corner, max_size);
738 void closest_pair(
unsigned int & ID1,
unsigned int & ID2,
739 double & distance2)
const;
740 void remove(
unsigned int ID);
742 virtual unsigned int replace(
unsigned int ID1,
unsigned int ID2,
744 virtual void replace_many(
const std::vector<unsigned int> & IDs_to_remove,
745 const std::vector<Coord2D> & new_positions,
746 std::vector<unsigned int> & new_IDs);
748 outdev <<
_trees[0]->max_depth() <<
" "
749 <<
_trees[1]->max_depth() <<
" "
750 <<
_trees[2]->max_depth() <<
"\n";
756 const unsigned int max_size);
809 if (x>y)
return false;
819 #endif // __FJCORE_CLOSESTPAIR2D__HH__
820 #ifndef __FJCORE_LAZYTILING9ALT_HH__
821 #define __FJCORE_LAZYTILING9ALT_HH__
868 return deta*deta + dphi*
dphi;
873 return deta*deta + dphi*
dphi;
878 return deta*deta + dphi*
dphi;
883 return deta*deta + dphi*
dphi;
892 const std::vector<PseudoJet> &
_jets;
910 std::vector<int> & tile_union,
int & n_near_tiles)
const;
912 std::vector<int> & tile_union,
int & n_near_tiles);
914 std::vector<int> & tile_union,
int & n_near_tiles);
916 void _set_NN(
TiledJet * jetI, std::vector<TiledJet *> & jets_for_minheap);
917 template <
class J>
double _bj_diJ(
const J *
const jet)
const {
918 double kt2 = jet->kt2;
919 if (jet->NN != NULL) {
if (jet->NN->kt2 < kt2) {kt2 = jet->NN->kt2;}}
920 return jet->NN_dist * kt2;
923 J *
const jetA,
const int _jets_index)
const {
924 jetA->eta =
_jets[_jets_index].rap();
925 jetA->phi =
_jets[_jets_index].phi_02pi();
926 jetA->kt2 =
_cs.jet_scale_for_algorithm(
_jets[_jets_index]);
927 jetA->_jets_index = _jets_index;
932 const J *
const jetA,
const J *
const jetB)
const {
933 double dphi = std::abs(jetA->phi - jetB->phi);
934 double deta = (jetA->eta - jetB->eta);
936 return dphi*dphi + deta*deta;
939 const J *
const jetA,
const J *
const jetB)
const {
940 double dphi = jetA->phi - jetB->phi;
941 double deta = (jetA->eta - jetB->eta);
942 return dphi*dphi + deta*deta;
946 #endif // __FJCORE_LAZYTILING9ALT_HH__
947 #ifndef __FJCORE_LAZYTILING9_HH__
948 #define __FJCORE_LAZYTILING9_HH__
979 const std::vector<PseudoJet> &
_jets;
984 #endif // INSTRUMENT2
1001 std::vector<int> & tile_union,
int & n_near_tiles)
const;
1003 std::vector<int> & tile_union,
int & n_near_tiles);
1005 std::vector<int> & tile_union,
int & n_near_tiles);
1013 void _set_NN(
TiledJet * jetI, std::vector<TiledJet *> & jets_for_minheap);
1014 template <
class J>
double _bj_diJ(
const J *
const jet)
const {
1015 double kt2 = jet->kt2;
1016 if (jet->NN != NULL) {
if (jet->NN->kt2 < kt2) {kt2 = jet->NN->kt2;}}
1017 return jet->NN_dist * kt2;
1020 J *
const jetA,
const int _jets_index)
const {
1021 jetA->eta =
_jets[_jets_index].rap();
1022 jetA->phi =
_jets[_jets_index].phi_02pi();
1023 jetA->kt2 =
_cs.jet_scale_for_algorithm(
_jets[_jets_index]);
1024 jetA->_jets_index = _jets_index;
1025 jetA->NN_dist =
_R2;
1029 const J *
const jetA,
const J *
const jetB)
1036 double dphi = std::abs(jetA->phi - jetB->phi);
1037 double deta = (jetA->eta - jetB->eta);
1039 return dphi*dphi + deta*deta;
1042 const J *
const jetA,
const J *
const jetB)
1049 double dphi = jetA->phi - jetB->phi;
1050 double deta = (jetA->eta - jetB->eta);
1051 return dphi*dphi + deta*deta;
1055 #endif // __FJCORE_LAZYTILING9_HH__
1056 #ifndef __FJCORE_LAZYTILING25_HH__
1057 #define __FJCORE_LAZYTILING25_HH__
1071 #endif // INSTRUMENT2
1088 std::vector<int> & tile_union,
int & n_near_tiles)
const;
1090 std::vector<int> & tile_union,
int & n_near_tiles);
1092 std::vector<int> & tile_union,
int & n_near_tiles);
1100 void _set_NN(
TiledJet * jetI, std::vector<TiledJet *> & jets_for_minheap);
1101 template <
class J>
double _bj_diJ(
const J *
const jet)
const {
1102 double kt2 = jet->kt2;
1103 if (jet->NN != NULL) {
if (jet->NN->kt2 < kt2) {kt2 = jet->NN->kt2;}}
1104 return jet->NN_dist * kt2;
1107 J *
const jetA,
const int _jets_index)
const {
1108 jetA->eta =
_jets[_jets_index].rap();
1109 jetA->phi =
_jets[_jets_index].phi_02pi();
1110 jetA->kt2 =
_cs.jet_scale_for_algorithm(
_jets[_jets_index]);
1111 jetA->_jets_index = _jets_index;
1112 jetA->NN_dist =
_R2;
1116 const J *
const jetA,
const J *
const jetB)
1123 double dphi = std::abs(jetA->phi - jetB->phi);
1124 double deta = (jetA->eta - jetB->eta);
1126 return dphi*dphi + deta*deta;
1129 const J *
const jetA,
const J *
const jetB)
1136 double dphi = jetA->phi - jetB->phi;
1137 double deta = (jetA->eta - jetB->eta);
1138 return dphi*dphi + deta*deta;
1142 #endif // __FJCORE_LAZYTILING25_HH__
1143 #ifndef __FJCORE_TILINGEXTENT_HH__
1144 #define __FJCORE_TILINGEXTENT_HH__
1155 void _determine_rapidity_extent(
const std::vector<PseudoJet> & particles);
1158 #endif // __FJCORE_TILINGEXTENT_HH__
1165 using namespace std;
1167 unsigned int shift) {
1168 Coord2D renorm_point = (point.
coord - _left_corner)/_range;
1173 shuffle.
x =
static_cast<unsigned int>(twopow31 * renorm_point.
x) + shift;
1174 shuffle.
y =
static_cast<unsigned int>(twopow31 * renorm_point.
y) + shift;
1187 unsigned int max_size) {
1188 unsigned int n_positions = positions.size();
1189 assert(max_size >= n_positions);
1190 _points.resize(max_size);
1191 for (
unsigned int i = n_positions;
i < max_size;
i++) {
1192 _available_points.push(&(_points[
i]));
1194 _left_corner = left_corner;
1195 _range = max((right_corner.
x - left_corner.
x),
1196 (right_corner.
y - left_corner.
y));
1197 vector<Shuffle> shuffles(n_positions);
1198 for (
unsigned int i = 0;
i < n_positions;
i++) {
1199 _points[
i].coord = positions[
i];
1200 _points[
i].neighbour_dist2 = numeric_limits<double>::max();
1201 _points[
i].review_flag = 0;
1202 _point2shuffle(_points[
i], shuffles[i], 0);
1204 for (
unsigned ishift = 0; ishift < _nshift; ishift++) {
1205 _shifts[ishift] =
static_cast<unsigned int>(((twopow31*1.0)*ishift)/_nshift);
1206 if (ishift == 0) {_rel_shifts[ishift] = 0;}
1207 else {_rel_shifts[ishift] = _shifts[ishift] - _shifts[ishift-1];}
1209 _cp_search_range = 30;
1210 _points_under_review.reserve(_nshift * _cp_search_range);
1211 for (
unsigned int ishift = 0; ishift < _nshift; ishift++) {
1213 unsigned rel_shift = _rel_shifts[ishift];
1214 for (
unsigned int i = 0;
i < shuffles.size();
i++) {
1215 shuffles[
i] += rel_shift; }
1217 sort(shuffles.begin(), shuffles.end());
1220 unsigned int CP_range =
min(_cp_search_range, n_positions-1);
1222 Point * this_point = circ->point;
1223 this_point->
circ[ishift] = circ;
1225 for (
unsigned i=0;
i < CP_range;
i++) {
1227 double dist2 = this_point->
distance2(*other->point);
1228 if (dist2 < this_point->neighbour_dist2) {
1233 }
while (++circ !=
start);
1235 vector<double> mindists2(n_positions);
1236 for (
unsigned int i = 0;
i < n_positions;
i++) {
1237 mindists2[
i] = _points[
i].neighbour_dist2;}
1241 double & distance2)
const {
1242 ID1 = _heap->minloc();
1243 ID2 = _ID(_points[ID1].neighbour);
1244 distance2 = _points[ID1].neighbour_dist2;
1248 if (point->
review_flag == 0) _points_under_review.push_back(point);
1252 if (point->
review_flag == 0) _points_under_review.push_back(point);
1256 Point * point_to_remove = & (_points[ID]);
1257 _remove_from_search_tree(point_to_remove);
1258 _deal_with_points_to_review();
1261 _available_points.push(point_to_remove);
1262 _set_label(point_to_remove, _remove_heap_entry);
1263 unsigned int CP_range =
min(_cp_search_range,
size()-1);
1264 for (
unsigned int ishift = 0; ishift < _nshift; ishift++) {
1267 _trees[ishift]->remove(removed_circ);
1268 circulator left_end = right_end, orig_right_end = right_end;
1269 for (
unsigned int i = 0;
i < CP_range;
i++) {left_end--;}
1270 if (
size()-1 < _cp_search_range) {
1271 left_end--; right_end--;
1274 Point * left_point = left_end->point;
1275 if (left_point->
neighbour == point_to_remove) {
1277 _add_label(left_point, _review_neighbour);
1280 double dist2 = left_point->
distance2(*right_end->point);
1281 if (dist2 < left_point->neighbour_dist2) {
1282 left_point->
neighbour = right_end->point;
1285 _add_label(left_point, _review_heap_entry);
1289 }
while (++left_end != orig_right_end);
1293 unsigned int CP_range =
min(_cp_search_range,
size()-1);
1294 while(_points_under_review.size() > 0) {
1295 Point * this_point = _points_under_review.back();
1296 _points_under_review.pop_back();
1297 if (this_point->
review_flag & _remove_heap_entry) {
1299 _heap->remove(_ID(this_point));
1302 if (this_point->
review_flag & _review_neighbour) {
1305 for (
unsigned int ishift = 0; ishift < _nshift; ishift++) {
1308 for (
unsigned i=0;
i < CP_range;
i++) {
1310 double dist2 = this_point->
distance2(*other->point);
1311 if (dist2 < this_point->neighbour_dist2) {
1324 assert(_available_points.size() > 0);
1325 Point * new_point = _available_points.top();
1326 _available_points.pop();
1327 new_point->
coord = new_coord;
1328 _insert_into_search_tree(new_point);
1329 _deal_with_points_to_review();
1330 return _ID(new_point);
1334 Point * point_to_remove = & (_points[ID1]);
1335 _remove_from_search_tree(point_to_remove);
1336 point_to_remove = & (_points[ID2]);
1337 _remove_from_search_tree(point_to_remove);
1338 Point * new_point = _available_points.top();
1339 _available_points.pop();
1341 _insert_into_search_tree(new_point);
1342 _deal_with_points_to_review();
1343 return _ID(new_point);
1346 const std::vector<unsigned int> & IDs_to_remove,
1347 const std::vector<Coord2D> & new_positions,
1348 std::vector<unsigned int> & new_IDs) {
1349 for (
unsigned int i = 0;
i < IDs_to_remove.size();
i++) {
1350 _remove_from_search_tree(& (_points[IDs_to_remove[
i]]));
1353 for (
unsigned int i = 0;
i < new_positions.size();
i++) {
1354 Point * new_point = _available_points.top();
1355 _available_points.pop();
1356 new_point->
coord = new_positions[
i];
1357 _insert_into_search_tree(new_point);
1358 new_IDs.push_back(_ID(new_point));
1360 _deal_with_points_to_review();
1363 _set_label(new_point, _review_heap_entry);
1365 unsigned int CP_range =
min(_cp_search_range,
size()-1);
1366 for (
unsigned ishift = 0; ishift < _nshift; ishift++) {
1368 _point2shuffle(*new_point, new_shuffle, _shifts[ishift]);
1369 circulator new_circ = _trees[ishift]->insert(new_shuffle);
1370 new_point->
circ[ishift] = new_circ;
1371 circulator right_edge = new_circ; right_edge++;
1373 for (
unsigned int i = 0;
i < CP_range;
i++) {left_edge--;}
1375 Point * left_point = left_edge->point;
1376 Point * right_point = right_edge->point;
1377 double new_dist2 = left_point->
distance2(*new_point);
1378 if (new_dist2 < left_point->neighbour_dist2) {
1381 _add_label(left_point, _review_heap_entry);
1383 new_dist2 = new_point->
distance2(*right_point);
1384 if (new_dist2 < new_point->neighbour_dist2) {
1388 if (left_point->
neighbour == right_point) {
1389 _add_label(left_point, _review_neighbour);
1392 }
while (++left_edge != new_circ);
1405 using namespace std;
1408 if (_structure_shared_ptr){
1412 if (_deletes_self_when_unused) {
1413 _structure_shared_ptr.set_count(_structure_shared_ptr.use_count()
1414 + _structure_use_count_after_construction);
1419 assert(_deletes_self_when_unused);
1420 _deletes_self_when_unused =
false;
1423 const JetDefinition & jet_def_in,
1424 const bool & writeout_combinations) {
1425 _decant_options(jet_def_in, writeout_combinations);
1426 _initialise_and_run_no_decant();
1429 _fill_initial_history();
1430 if (n_particles() == 0)
return;
1432 _plugin_activated =
true;
1433 _jet_def.plugin()->run_clustering( (*
this) );
1434 _plugin_activated =
false;
1435 _update_structure_use_count();
1454 _simple_N2_cluster_EEBriefJet();
1457 throw Error(
"A ClusterSequence cannot be created with an uninitialised JetDefinition");
1459 if (_strategy ==
Best) {
1460 _strategy = _best_strategy();
1461 #ifdef __FJCORE_DROP_CGAL
1463 #endif // __FJCORE_DROP_CGAL
1464 }
else if (_strategy ==
BestFJ30) {
1470 #ifndef __FJCORE_DROP_CGAL
1472 || N > 35000/pow(
_Rparam,1.15)) {
1474 #endif // __FJCORE_DROP_CGAL
1475 }
else if (N <= 450) {
1482 if ( _strategy ==
NlnN
1487 #ifdef __FJCORE_DROP_CGAL
1493 if (_jet_def.strategy() !=
Best && _strategy != _jet_def.strategy()) {
1495 oss <<
"Cluster strategy " << strategy_string(_jet_def.strategy())
1496 <<
" automatically changed to " << strategy_string()
1497 <<
" because the former is not supported for R = " <<
_Rparam
1499 _changed_strategy_warning.warn(oss.str());
1503 this->_simple_N2_cluster_BriefJet();
1504 }
else if (_strategy ==
N2Tiled) {
1505 this->_faster_tiled_N2_cluster();
1507 this->_minheap_faster_tiled_N2_cluster();
1509 _plugin_activated =
true;
1512 _plugin_activated =
false;
1514 _plugin_activated =
true;
1517 _plugin_activated =
false;
1519 _plugin_activated =
true;
1522 _plugin_activated =
false;
1524 throw Error(
"N2MHTLazy9AntiKtSeparateGhosts strategy not supported with FJCORE");
1525 }
else if (_strategy ==
NlnN) {
1526 this->_delaunay_cluster();
1527 }
else if (_strategy ==
NlnNCam) {
1528 this->_CP2DChan_cluster_2piMultD();
1530 this->_delaunay_cluster();
1531 }
else if (_strategy ==
N3Dumb ) {
1532 this->_really_dumb_cluster();
1534 this->_tiled_N2_cluster();
1536 this->_CP2DChan_cluster();
1538 this->_CP2DChan_cluster_2pi2R();
1541 err <<
"Unrecognised value for strategy: "<<_strategy;
1542 throw Error(err.str());
1551 if (!_first_time) {
return;}
1552 _first_time =
false;
1553 ostream * ostr = _fastjet_banner_ostr;
1555 (*ostr) <<
"#--------------------------------------------------------------------------\n";
1556 (*ostr) <<
"# FastJet release " <<
fastjet_version <<
" [fjcore]" << endl;
1557 (*ostr) <<
"# M. Cacciari, G.P. Salam and G. Soyez \n";
1558 (*ostr) <<
"# A software package for jet finding and analysis at colliders \n";
1559 (*ostr) <<
"# http://fastjet.fr \n";
1561 (*ostr) <<
"# Please cite EPJC72(2012)1896 [arXiv:1111.6097] if you use this package\n";
1562 (*ostr) <<
"# for scientific work and optionally PLB641(2006)57 [hep-ph/0512210]. \n";
1564 (*ostr) <<
"# FastJet is provided without warranty under the terms of the GNU GPLv2.\n";
1565 (*ostr) <<
"# It uses T. Chan's closest pair algorithm, S. Fortune's Voronoi code";
1566 #ifndef __FJCORE_DROP_CGAL
1567 (*ostr) <<
",\n# CGAL ";
1570 #endif // __FJCORE_DROP_CGAL
1571 (*ostr) <<
"and 3rd party plugin jet algorithms. See COPYING file for details.\n";
1572 (*ostr) <<
"#--------------------------------------------------------------------------\n";
1576 const bool & writeout_combinations) {
1577 _jet_def = jet_def_in;
1578 _writeout_combinations = writeout_combinations;
1580 _decant_options_partial();
1584 _jet_algorithm = _jet_def.jet_algorithm();
1586 _strategy = _jet_def.strategy();
1587 _plugin_activated =
false;
1588 _update_structure_use_count();
1592 _history.reserve(
_jets.size()*2);
1594 for (
int i = 0; i < static_cast<int>(
_jets.size()) ;
i++) {
1596 element.
parent1 = InexistentParent;
1597 element.
parent2 = InexistentParent;
1602 _history.push_back(element);
1603 _jet_def.recombiner()->preprocess(
_jets[
i]);
1604 _jets[
i].set_cluster_hist_index(i);
1605 _set_structure_shared_ptr(
_jets[i]);
1608 _initial_n =
_jets.size();
1609 _deletes_self_when_unused =
false;
1613 switch(strategy_in) {
1615 strategy =
"NlnN";
break;
1617 strategy =
"NlnN3pi";
break;
1619 strategy =
"NlnN4pi";
break;
1621 strategy =
"N2Plain";
break;
1623 strategy =
"N2Tiled";
break;
1625 strategy =
"N2MinHeapTiled";
break;
1627 strategy =
"N2PoorTiled";
break;
1629 strategy =
"N2MHTLazy9";
break;
1631 strategy =
"N2MHTLazy9Alt";
break;
1633 strategy =
"N2MHTLazy25";
break;
1635 strategy =
"N2MHTLazy9AntiKtSeparateGhosts";
break;
1637 strategy =
"N3Dumb";
break;
1639 strategy =
"NlnNCam4pi";
break;
1641 strategy =
"NlnNCam2pi2R";
break;
1643 strategy =
"NlnNCam";
break;
1645 strategy =
"plugin strategy";
break;
1647 strategy =
"Unrecognized";
1652 const PseudoJet & jet)
const {
1653 if (_jet_algorithm ==
kt_algorithm) {
return jet.kt2();}
1656 double kt2=jet.kt2();
1657 return kt2 > 1
e-300 ? 1.0/kt2 : 1e300;
1659 double kt2 = jet.kt2();
1660 double p = jet_def().extra_param();
1661 if (p <= 0 && kt2 < 1
e-300) kt2 = 1
e-300;
1664 double kt2 = jet.kt2();
1665 double lim = _jet_def.extra_param();
1666 if (kt2 < lim*lim && kt2 != 0.0) {
1668 }
else {
return 1.0;}
1669 }
else {
throw Error(
"Unrecognised jet algorithm");}
1673 double bounded_R = max(
_Rparam, 0.1);
1674 if (N <= 30 || N <= 39.0/(bounded_R + 0.6)) {
1677 const static _Parabola N_Tiled_to_MHT_lowR (-45.4947,54.3528,44.6283);
1678 const static _Parabola L_MHT_to_MHTLazy9_lowR (0.677807,-1.05006,10.6994);
1679 const static _Parabola L_MHTLazy9_to_MHTLazy25_akt_lowR(0.169967,-0.512589,12.1572);
1680 const static _Parabola L_MHTLazy9_to_MHTLazy25_kt_lowR (0.16237,-0.484612,12.3373);
1681 const static _Parabola L_MHTLazy9_to_MHTLazy25_cam_lowR = L_MHTLazy9_to_MHTLazy25_kt_lowR;
1682 const static _Parabola L_MHTLazy25_to_NlnN_akt_lowR (0.0472051,-0.22043,15.9196);
1683 const static _Parabola L_MHTLazy25_to_NlnN_kt_lowR (0.118609,-0.326811,14.8287);
1684 const static _Parabola L_MHTLazy25_to_NlnN_cam_lowR (0.10119,-0.295748,14.3924);
1685 const static _Line L_Tiled_to_MHTLazy9_medR (-1.31304,7.29621);
1686 const static _Parabola L_MHTLazy9_to_MHTLazy25_akt_medR = L_MHTLazy9_to_MHTLazy25_akt_lowR;
1687 const static _Parabola L_MHTLazy9_to_MHTLazy25_kt_medR = L_MHTLazy9_to_MHTLazy25_kt_lowR;
1688 const static _Parabola L_MHTLazy9_to_MHTLazy25_cam_medR = L_MHTLazy9_to_MHTLazy25_cam_lowR;
1689 const static _Parabola L_MHTLazy25_to_NlnN_akt_medR = L_MHTLazy25_to_NlnN_akt_lowR;
1690 const static _Parabola L_MHTLazy25_to_NlnN_kt_medR = L_MHTLazy25_to_NlnN_kt_lowR;
1691 const static _Parabola L_MHTLazy25_to_NlnN_cam_medR = L_MHTLazy25_to_NlnN_cam_lowR;
1692 const static double N_Plain_to_MHTLazy9_largeR = 75;
1693 const static double N_MHTLazy9_to_MHTLazy25_akt_largeR = 700;
1694 const static double N_MHTLazy9_to_MHTLazy25_kt_largeR = 1000;
1695 const static double N_MHTLazy9_to_MHTLazy25_cam_largeR = 1000;
1696 const static double N_MHTLazy25_to_NlnN_akt_largeR = 100000;
1697 const static double N_MHTLazy25_to_NlnN_kt_largeR = 40000;
1698 const static double N_MHTLazy25_to_NlnN_cam_largeR = 15000;
1701 double p = jet_def().extra_param();
1707 jet_algorithm = _jet_algorithm;
1709 if (bounded_R < 0.65) {
1710 if (N < N_Tiled_to_MHT_lowR(bounded_R))
return N2Tiled;
1711 double logN = log(
double(N));
1712 if (logN < L_MHT_to_MHTLazy9_lowR(bounded_R))
return N2MinHeapTiled;
1715 if (logN < L_MHTLazy9_to_MHTLazy25_akt_lowR(bounded_R))
return N2MHTLazy9;
1716 else if (logN < L_MHTLazy25_to_NlnN_akt_lowR(bounded_R))
return N2MHTLazy25;
1719 if (logN < L_MHTLazy9_to_MHTLazy25_kt_lowR(bounded_R))
return N2MHTLazy9;
1720 else if (logN < L_MHTLazy25_to_NlnN_kt_lowR(bounded_R))
return N2MHTLazy25;
1723 if (logN < L_MHTLazy9_to_MHTLazy25_cam_lowR(bounded_R))
return N2MHTLazy9;
1724 else if (logN < L_MHTLazy25_to_NlnN_cam_lowR(bounded_R))
return N2MHTLazy25;
1728 }
else if (bounded_R < 0.5*
pi) {
1729 double logN = log(
double(N));
1730 if (logN < L_Tiled_to_MHTLazy9_medR(bounded_R))
return N2Tiled;
1733 if (logN < L_MHTLazy9_to_MHTLazy25_akt_medR(bounded_R))
return N2MHTLazy9;
1734 else if (logN < L_MHTLazy25_to_NlnN_akt_medR(bounded_R))
return N2MHTLazy25;
1737 if (logN < L_MHTLazy9_to_MHTLazy25_kt_medR(bounded_R))
return N2MHTLazy9;
1738 else if (logN < L_MHTLazy25_to_NlnN_kt_medR(bounded_R))
return N2MHTLazy25;
1741 if (logN < L_MHTLazy9_to_MHTLazy25_cam_medR(bounded_R))
return N2MHTLazy9;
1742 else if (logN < L_MHTLazy25_to_NlnN_cam_medR(bounded_R))
return N2MHTLazy25;
1747 if (N < N_Plain_to_MHTLazy9_largeR)
return N2Plain;
1750 if (N < N_MHTLazy9_to_MHTLazy25_akt_largeR)
return N2MHTLazy9;
1751 else if (N < N_MHTLazy25_to_NlnN_akt_largeR)
return N2MHTLazy25;
1754 if (N < N_MHTLazy9_to_MHTLazy25_kt_largeR)
return N2MHTLazy9;
1755 else if (N < N_MHTLazy25_to_NlnN_kt_largeR)
return N2MHTLazy25;
1758 if (N < N_MHTLazy9_to_MHTLazy25_cam_largeR)
return N2MHTLazy9;
1759 else if (N < N_MHTLazy25_to_NlnN_cam_largeR)
return N2MHTLazy25;
1764 assert(0 &&
"Code should never reach here");
1769 _deletes_self_when_unused =
false;
1770 transfer_from_sequence(cs);
1776 if (will_delete_self_when_unused())
1777 throw(
Error(
"cannot use CS::transfer_from_sequence after a call to delete_self_when_unused()"));
1778 _jet_def = from_seq._jet_def ;
1779 _writeout_combinations = from_seq._writeout_combinations ;
1780 _initial_n = from_seq._initial_n ;
1782 _R2 = from_seq._R2 ;
1783 _invR2 = from_seq._invR2 ;
1784 _strategy = from_seq._strategy ;
1785 _jet_algorithm = from_seq._jet_algorithm ;
1786 _plugin_activated = from_seq._plugin_activated ;
1788 _jets = (*action_on_jets)(from_seq._jets);
1790 _jets = from_seq._jets;
1791 _history = from_seq._history;
1792 _extras = from_seq._extras;
1793 if (_structure_shared_ptr) {
1794 if (_deletes_self_when_unused)
throw Error(
"transfer_from_sequence cannot be used for a cluster sequence that deletes self when unused");
1800 _update_structure_use_count();
1801 for (
unsigned int i=0;
i<_jets.size();
i++){
1802 _jets[
i].set_cluster_hist_index(from_seq._jets[
i].cluster_hist_index());
1803 _set_structure_shared_ptr(_jets[
i]);
1807 int jet_i,
int jet_j,
double dij,
1808 const PseudoJet & newjet,
int & newjet_k) {
1809 plugin_record_ij_recombination(jet_i, jet_j, dij, newjet_k);
1810 int tmp_index =
_jets[newjet_k].cluster_hist_index();
1811 _jets[newjet_k] = newjet;
1812 _jets[newjet_k].set_cluster_hist_index(tmp_index);
1813 _set_structure_shared_ptr(
_jets[newjet_k]);
1816 double dcut = ptmin*ptmin;
1817 int i = _history.size() - 1;
1818 vector<PseudoJet> jets_local;
1821 if (_history[i].max_dij_so_far < dcut) {
break;}
1822 if (_history[i].parent2 == BeamJet && _history[i].dij >= dcut) {
1824 int parent1 = _history[
i].parent1;
1825 jets_local.push_back(
_jets[_history[parent1].jetp_index]);}
1830 if (_history[i].parent2 != BeamJet) {
break;}
1831 int parent1 = _history[
i].parent1;
1832 const PseudoJet & jet =
_jets[_history[parent1].jetp_index];
1833 if (jet.perp2() >= dcut) {jets_local.push_back(jet);}
1843 if (_history[i].parent2 == BeamJet) {
1844 int parent1 = _history[
i].parent1;
1845 const PseudoJet & jet =
_jets[_history[parent1].jetp_index];
1846 if (jet.perp2() >= dcut) {jets_local.push_back(jet);}
1850 }
else {
throw Error(
"cs::inclusive_jets(...): Unrecognized jet algorithm");}
1854 int i = _history.size() - 1;
1856 if (_history[i].max_dij_so_far <= dcut) {
break;}
1859 int stop_point = i + 1;
1860 int njets = 2*_initial_n - stop_point;
1864 int njets = n_exclusive_jets(dcut);
1865 return exclusive_jets(njets);
1868 if (njets > _initial_n) {
1870 err <<
"Requested " << njets <<
" exclusive jets, but there were only "
1871 << _initial_n <<
" particles in the event";
1872 throw Error(err.str());
1874 return exclusive_jets_up_to(njets);
1882 (_jet_def.extra_param() <0)) &&
1884 (!_jet_def.plugin()->exclusive_sequence_meaningful()))) {
1885 _exclusive_warnings.warn(
"dcut and exclusive jets for jet-finders other than kt, C/A or genkt with p>=0 should be interpreted with care.");
1887 int stop_point = 2*_initial_n - njets;
1888 if (stop_point < _initial_n) stop_point = _initial_n;
1889 if (2*_initial_n != static_cast<int>(_history.size())) {
1891 err <<
"2*_initial_n != _history.size() -- this endangers internal assumptions!\n";
1892 throw Error(err.str());
1894 vector<PseudoJet> jets_local;
1895 for (
unsigned int i = stop_point;
i < _history.size();
i++) {
1896 int parent1 = _history[
i].parent1;
1897 if (parent1 < stop_point) {
1898 jets_local.push_back(
_jets[_history[parent1].jetp_index]);
1900 int parent2 = _history[
i].parent2;
1901 if (parent2 < stop_point && parent2 > 0) {
1902 jets_local.push_back(
_jets[_history[parent2].jetp_index]);
1905 if (
int(jets_local.size()) !=
min(_initial_n, njets)) {
1907 err <<
"ClusterSequence::exclusive_jets: size of returned vector ("
1908 <<jets_local.size()<<
") does not coincide with requested number of jets ("
1910 throw Error(err.str());
1916 if (njets >= _initial_n) {
return 0.0;}
1917 return _history[2*_initial_n-njets-1].dij;
1921 if (njets >= _initial_n) {
return 0.0;}
1922 return _history[2*_initial_n-njets-1].max_dij_so_far;
1925 (
const PseudoJet & jet,
const double dcut)
const {
1926 set<const history_element*> subhist;
1927 get_subhist_set(subhist, jet, dcut, 0);
1928 vector<PseudoJet> subjets;
1929 subjets.reserve(subhist.size());
1930 for (set<const history_element*>::iterator elem = subhist.begin();
1931 elem != subhist.end(); elem++) {
1932 subjets.push_back(
_jets[(*elem)->jetp_index]);
1937 const double dcut)
const {
1938 set<const history_element*> subhist;
1939 get_subhist_set(subhist, jet, dcut, 0);
1940 return subhist.size();
1943 (
const PseudoJet & jet,
int nsub)
const {
1944 vector<PseudoJet> subjets = exclusive_subjets_up_to(jet, nsub);
1945 if (
int(subjets.size()) < nsub) {
1947 err <<
"Requested " << nsub <<
" exclusive subjets, but there were only "
1948 << subjets.size() <<
" particles in the jet";
1949 throw Error(err.str());
1954 (
const PseudoJet & jet,
int nsub)
const {
1955 set<const history_element*> subhist;
1956 vector<PseudoJet> subjets;
1957 if (nsub < 0)
throw Error(
"Requested a negative number of subjets. This is nonsensical.");
1958 if (nsub == 0)
return subjets;
1959 get_subhist_set(subhist, jet, -1.0, nsub);
1960 subjets.reserve(subhist.size());
1961 for (set<const history_element*>::iterator elem = subhist.begin();
1962 elem != subhist.end(); elem++) {
1963 subjets.push_back(
_jets[(*elem)->jetp_index]);
1968 set<const history_element*> subhist;
1969 get_subhist_set(subhist, jet, -1.0, nsub);
1970 set<const history_element*>::iterator highest = subhist.end();
1972 return (*highest)->dij;
1975 set<const history_element*> subhist;
1976 get_subhist_set(subhist, jet, -1.0, nsub);
1977 set<const history_element*>::iterator highest = subhist.end();
1979 return (*highest)->max_dij_so_far;
1982 const PseudoJet & jet,
1983 double dcut,
int maxjet)
const {
1986 subhist.insert(&(_history[jet.cluster_hist_index()]));
1989 set<const history_element*>::iterator highest = subhist.end();
1990 assert (highest != subhist.begin());
1993 if (njet == maxjet)
break;
1996 subhist.erase(highest);
1997 subhist.insert(&(_history[elem->
parent1]));
1998 subhist.insert(&(_history[elem->
parent2]));
2003 const PseudoJet & jet)
const {
2004 assert(contains(
object) && contains(jet));
2005 const PseudoJet * this_object = &object;
2006 const PseudoJet * childp;
2008 if (this_object->cluster_hist_index() == jet.cluster_hist_index()) {
2010 }
else if (has_child(*this_object, childp)) {
2011 this_object = childp;
2018 PseudoJet & parent2)
const {
2023 parent1 = PseudoJet(0.0,0.0,0.0,0.0);
2029 if (parent1.perp2() < parent2.perp2())
std::swap(parent1,parent2);
2034 const PseudoJet * childp;
2035 bool res = has_child(jet, childp);
2040 child = PseudoJet(0.0,0.0,0.0,0.0);
2045 const history_element &
hist = _history[jet.cluster_hist_index()];
2046 if (hist.child >= 0 && _history[hist.child].jetp_index >= 0) {
2047 childp = &(
_jets[_history[hist.child].jetp_index]);
2055 PseudoJet & partner)
const {
2057 if (hist.
child >= 0 && _history[hist.
child].parent2 >= 0) {
2059 if (child_hist.
parent1 == jet.cluster_hist_index()) {
2060 partner =
_jets[_history[child_hist.
parent2].jetp_index];
2062 partner =
_jets[_history[child_hist.
parent1].jetp_index];
2066 partner = PseudoJet(0.0,0.0,0.0,0.0);
2071 vector<PseudoJet> subjets;
2072 add_constituents(jet, subjets);
2076 ostream & ostr)
const {
2077 for (
unsigned i = 0;
i < jets_in.size();
i++) {
2079 << jets_in[
i].px() <<
" "
2080 << jets_in[
i].py() <<
" "
2081 << jets_in[
i].pz() <<
" "
2082 << jets_in[
i].E() << endl;
2083 vector<PseudoJet> cst = constituents(jets_in[
i]);
2084 for (
unsigned j = 0;
j < cst.size() ;
j++) {
2085 ostr <<
" " <<
j <<
" "
2086 << cst[
j].rap() <<
" "
2087 << cst[
j].phi() <<
" "
2088 << cst[
j].perp() << endl;
2090 ostr <<
"#END" << endl;
2096 std::ofstream ostr(filename.c_str());
2097 if (comment !=
"") ostr <<
"# " << comment << endl;
2098 print_jets_for_root(jets_in, ostr);
2101 const vector<PseudoJet> & jets_in)
const {
2102 vector<int> indices(n_particles());
2103 for (
unsigned ipart = 0; ipart < n_particles(); ipart++)
2104 indices[ipart] = -1;
2105 for (
unsigned ijet = 0; ijet < jets_in.size(); ijet++) {
2106 vector<PseudoJet> jet_constituents(constituents(jets_in[ijet]));
2107 for (
unsigned ip = 0; ip < jet_constituents.size(); ip++) {
2108 unsigned iclust = jet_constituents[ip].cluster_hist_index();
2109 unsigned ipart = history()[iclust].jetp_index;
2110 indices[ipart] = ijet;
2116 const PseudoJet & jet, vector<PseudoJet> & subjet_vector)
const {
2117 int i = jet.cluster_hist_index();
2118 int parent1 = _history[
i].parent1;
2119 int parent2 = _history[
i].parent2;
2120 if (parent1 == InexistentParent) {
2121 subjet_vector.push_back(
_jets[i]);
2124 add_constituents(
_jets[_history[parent1].jetp_index], subjet_vector);
2125 if (parent2 != BeamJet) {
2126 add_constituents(
_jets[_history[parent2].jetp_index], subjet_vector);
2131 const int parent2,
const int jetp_index,
2139 element.
max_dij_so_far = max(dij,_history[_history.size()-1].max_dij_so_far);
2140 _history.push_back(element);
2141 int local_step = _history.size()-1;
2143 if (_history[parent1].child !=
Invalid){
2144 throw InternalError(
"trying to recomine an object that has previsously been recombined");
2146 _history[parent1].child = local_step;
2148 if (_history[parent2].child !=
Invalid){
2149 throw InternalError(
"trying to recomine an object that has previsously been recombined");
2151 _history[parent2].child = local_step;
2155 _jets[jetp_index].set_cluster_hist_index(local_step);
2156 _set_structure_shared_ptr(
_jets[jetp_index]);
2158 if (_writeout_combinations) {
2159 cout << local_step <<
": "
2160 << parent1 <<
" with " << parent2
2161 <<
"; y = "<< dij<<endl;
2165 valarray<int> lowest_constituent(_history.size());
2166 int hist_n = _history.size();
2167 lowest_constituent = hist_n;
2168 for (
int i = 0;
i < hist_n;
i++) {
2169 lowest_constituent[
i] =
min(lowest_constituent[
i],i);
2170 if (_history[i].child > 0) lowest_constituent[_history[
i].child]
2171 =
min(lowest_constituent[_history[i].child],lowest_constituent[i]);
2173 valarray<bool> extracted(_history.size()); extracted =
false;
2174 vector<int> unique_tree;
2175 unique_tree.reserve(_history.size());
2176 for (
unsigned i = 0;
i < n_particles();
i++) {
2177 if (!extracted[
i]) {
2178 unique_tree.push_back(i);
2179 extracted[
i] =
true;
2180 _extract_tree_children(i, extracted, lowest_constituent, unique_tree);
2187 valarray<bool> & extracted,
2188 const valarray<int> & lowest_constituent,
2189 vector<int> & unique_tree)
const {
2190 if (!extracted[position]) {
2191 _extract_tree_parents(position,extracted,lowest_constituent,unique_tree);
2193 int child = _history[
position].child;
2194 if (child >= 0) _extract_tree_children(child,extracted,lowest_constituent,unique_tree);
2197 vector<PseudoJet> unclustered;
2198 for (
unsigned i = 0;
i < n_particles() ;
i++) {
2200 unclustered.push_back(
_jets[_history[
i].jetp_index]);
2205 vector<PseudoJet> unclustered;
2206 for (
unsigned i = 0;
i < _history.size() ;
i++) {
2207 if ((_history[
i].child ==
Invalid) && (_history[
i].parent2 != BeamJet))
2208 unclustered.push_back(
_jets[_history[
i].jetp_index]);
2213 return jet.cluster_hist_index() >= 0
2214 && jet.cluster_hist_index() < int(_history.size())
2215 && jet.has_valid_cluster_sequence()
2216 && jet.associated_cluster_sequence() ==
this;
2220 valarray<bool> & extracted,
2221 const valarray<int> & lowest_constituent,
2222 vector<int> & unique_tree)
const {
2223 if (!extracted[position]) {
2224 int parent1 = _history[
position].parent1;
2225 int parent2 = _history[
position].parent2;
2226 if (parent1 >= 0 && parent2 >= 0) {
2227 if (lowest_constituent[parent1] > lowest_constituent[parent2])
2230 if (parent1 >= 0 && !extracted[parent1])
2231 _extract_tree_parents(parent1,extracted,lowest_constituent,unique_tree);
2232 if (parent2 >= 0 && !extracted[parent2])
2233 _extract_tree_parents(parent2,extracted,lowest_constituent,unique_tree);
2234 unique_tree.push_back(position);
2239 const int jet_i,
const int jet_j,
2242 PseudoJet newjet(
false);
2243 _jet_def.recombiner()->recombine(
_jets[jet_i],
_jets[jet_j], newjet);
2244 _jets.push_back(newjet);
2245 newjet_k =
_jets.size()-1;
2246 int newstep_k = _history.size();
2247 _jets[newjet_k].set_cluster_hist_index(newstep_k);
2248 int hist_i =
_jets[jet_i].cluster_hist_index();
2249 int hist_j =
_jets[jet_j].cluster_hist_index();
2250 _add_step_to_history(
min(hist_i, hist_j), max(hist_i,hist_j),
2254 const int jet_i,
const double diB) {
2255 _add_step_to_history(
_jets[jet_i].cluster_hist_index(),BeamJet,
2260 j.set_structure_shared_ptr(_structure_shared_ptr);
2261 _update_structure_use_count();
2264 _structure_use_count_after_construction = _structure_shared_ptr.use_count();
2267 int new_count = _structure_shared_ptr.use_count() - _structure_use_count_after_construction;
2268 if (new_count <= 0) {
2269 throw Error(
"delete_self_when_unused may only be called if at least one object outside the CS (e.g. a jet) is already associated with the CS");
2271 _structure_shared_ptr.set_count(new_count);
2272 _deletes_self_when_unused =
true;
2280 using namespace std;
2289 if (point.
y < Dlim) {point.
y +=
twopi;
return true;}
2290 if (
twopi-point.
y < Dlim) {point.
y -=
twopi;
return true;}
2294 using namespace Private;
2296 unsigned int n = _initial_n;
2297 vector<MirrorInfo> coordIDs(2*n);
2298 vector<int> jetIDs(2*n);
2299 vector<Coord2D> coords(2*n);
2300 double Dlim4mirror =
min(Dlim,
pi);
2301 double minrap = numeric_limits<double>::max();
2302 double maxrap = -minrap;
2303 int coord_index = -1;
2305 for (
unsigned jet_i = 0; jet_i <
_jets.size(); jet_i++) {
2306 if (_history[
_jets[jet_i].cluster_hist_index()].child !=
Invalid ||
2307 (
_jets[jet_i].E() == abs(
_jets[jet_i].pz()) &&
2308 _jets[jet_i].perp2() == 0.0)
2311 coordIDs[jet_i].orig = ++coord_index;
2313 jetIDs[coord_index] = jet_i;
2314 minrap =
min(coords[coord_index].
x,minrap);
2315 maxrap = max(coords[coord_index].x,maxrap);
2316 Coord2D mirror_point(coords[coord_index]);
2318 coordIDs[jet_i].mirror = ++coord_index;
2319 coords[coord_index] = mirror_point;
2320 jetIDs[coord_index] = jet_i;
2322 coordIDs[jet_i].mirror =
Invalid;
2325 coords.resize(coord_index+1);
2326 Coord2D left_edge(minrap-1.0, -3.15);
2327 Coord2D right_edge(maxrap+1.0, 9.45);
2329 vector<Coord2D> new_points(2);
2330 vector<unsigned int> cIDs_to_remove(4);
2331 vector<unsigned int> new_cIDs(2);
2333 unsigned int cID1, cID2;
2336 if (distance2 > Dlim*Dlim) {
break;}
2338 int jet_i = jetIDs[cID1];
2339 int jet_j = jetIDs[cID2];
2342 _do_ij_recombination_step(jet_i, jet_j, distance2, newjet_k);
2343 if (--n_active == 1) {
break;}
2344 cIDs_to_remove.resize(0);
2345 cIDs_to_remove.push_back(coordIDs[jet_i].orig);
2346 cIDs_to_remove.push_back(coordIDs[jet_j].orig);
2347 if (coordIDs[jet_i].mirror !=
Invalid)
2348 cIDs_to_remove.push_back(coordIDs[jet_i].mirror);
2349 if (coordIDs[jet_j].mirror !=
Invalid)
2350 cIDs_to_remove.push_back(coordIDs[jet_j].mirror);
2352 new_points.resize(0);
2353 new_points.push_back(new_point);
2354 if (
make_mirror(new_point, Dlim4mirror)) new_points.push_back(new_point);
2356 coordIDs[newjet_k].orig = new_cIDs[0];
2357 jetIDs[new_cIDs[0]] = newjet_k;
2358 if (new_cIDs.size() == 2) {
2359 coordIDs[newjet_k].mirror = new_cIDs[1];
2360 jetIDs[new_cIDs[1]] = newjet_k;
2361 }
else {coordIDs[newjet_k].mirror =
Invalid;}
2365 if (_jet_algorithm !=
cambridge_algorithm)
throw Error(
"CP2DChan clustering method called for a jet-finder that is not the cambridge algorithm");
2366 _CP2DChan_limited_cluster(
_Rparam);
2367 _do_Cambridge_inclusive_jets();
2371 _CP2DChan_limited_cluster(
min(
_Rparam/2,0.3));
2373 _CP2DChan_cluster_2pi2R ();
2376 if (_jet_algorithm !=
cambridge_algorithm)
throw Error(
"_CP2DChan_cluster called for a jet-finder that is not the cambridge algorithm");
2377 unsigned int n =
_jets.size();
2378 vector<MirrorInfo> coordIDs(2*n);
2379 vector<int> jetIDs(2*n);
2380 vector<Coord2D> coords(2*n);
2381 double minrap = numeric_limits<double>::max();
2382 double maxrap = -minrap;
2383 int coord_index = 0;
2384 for (
unsigned i = 0;
i <
n;
i++) {
2388 coordIDs[
i].orig = coord_index;
2389 coordIDs[
i].mirror = coord_index+1;
2392 jetIDs[coord_index] =
i;
2393 jetIDs[coord_index+1] =
i;
2394 minrap =
min(coords[coord_index].
x,minrap);
2395 maxrap = max(coords[coord_index].x,maxrap);
2399 for (
unsigned i = n;
i < 2*
n;
i++) {coordIDs[
i].orig =
Invalid;}
2400 coords.resize(coord_index);
2401 Coord2D left_edge(minrap-1.0, 0.0);
2404 vector<Coord2D> new_points(2);
2405 vector<unsigned int> cIDs_to_remove(4);
2406 vector<unsigned int> new_cIDs(2);
2408 unsigned int cID1, cID2;
2412 if (distance2 > 1.0) {
break;}
2413 int jet_i = jetIDs[cID1];
2414 int jet_j = jetIDs[cID2];
2417 _do_ij_recombination_step(jet_i, jet_j, distance2, newjet_k);
2418 cIDs_to_remove[0] = coordIDs[jet_i].orig;
2419 cIDs_to_remove[1] = coordIDs[jet_i].mirror;
2420 cIDs_to_remove[2] = coordIDs[jet_j].orig;
2421 cIDs_to_remove[3] = coordIDs[jet_j].mirror;
2424 new_cIDs[0] = cp.
replace(cIDs_to_remove[0], cIDs_to_remove[2], new_points[0]);
2425 new_cIDs[1] = cp.
replace(cIDs_to_remove[1], cIDs_to_remove[3], new_points[1]);
2426 coordIDs[jet_i].orig =
Invalid;
2427 coordIDs[jet_j].orig =
Invalid;
2428 coordIDs[newjet_k] =
MirrorInfo(new_cIDs[0], new_cIDs[1]);
2429 jetIDs[new_cIDs[0]] = newjet_k;
2430 jetIDs[new_cIDs[1]] = newjet_k;
2432 if (n == 1) {
break;}
2434 _do_Cambridge_inclusive_jets();
2437 unsigned int n = _history.size();
2438 for (
unsigned int hist_i = 0; hist_i <
n; hist_i++) {
2439 if (_history[hist_i].child ==
Invalid) {
2440 _do_iB_recombination_step(_history[hist_i].jetp_index, 1.0);
2451 #ifndef __FJCORE_DROP_CGAL // in case we do not have the code for CGAL
2452 #include "fastjet/internal/Dnn4piCylinder.hh"
2453 #include "fastjet/internal/Dnn3piCylinder.hh"
2454 #include "fastjet/internal/Dnn2piCylinder.hh"
2455 #endif // __FJCORE_DROP_CGAL
2457 using namespace std;
2460 vector<EtaPhi> points(n);
2461 for (
int i = 0;
i <
n;
i++) {
2463 points[
i].sanitize();
2467 #ifndef __FJCORE_DROP_CGAL // strategy = NlnN* are not supported if we drop CGAL...
2470 DNN.
reset(
new Dnn4piCylinder(points,verbose));
2471 }
else if (_strategy ==
NlnN3pi) {
2472 DNN.
reset(
new Dnn3piCylinder(points,ignore_nearest_is_mirror,verbose));
2473 }
else if (_strategy ==
NlnN) {
2474 DNN.
reset(
new Dnn2piCylinder(points,ignore_nearest_is_mirror,verbose));
2479 err <<
"ERROR: Requested strategy "<<strategy_string()<<
" but it is not"<<endl;
2480 err <<
" supported because FastJet was compiled without CGAL"<<endl;
2481 throw Error(err.str());
2483 #endif // __FJCORE_DROP_CGAL
2488 for (
int ii = 0; ii <
n; ii++) {
2489 _add_ktdistance_to_map(ii, DijMap, DNN.
get());
2491 for (
int i=0;
i<
n;
i++) {
2496 bool recombine_with_beam;
2498 SmallestDij = DijMap.begin()->first;
2499 SmallestDijPair = DijMap.begin()->second;
2500 jet_i = SmallestDijPair.first;
2501 jet_j = SmallestDijPair.second;
2502 if (verbose) cout <<
"CS_Delaunay found recombination candidate: " << jet_i <<
" " << jet_j <<
" " << SmallestDij << endl;
2503 DijMap.erase(DijMap.begin());
2504 recombine_with_beam = (jet_j == BeamJet);
2505 if (!recombine_with_beam) {Valid2 = DNN->
Valid(jet_j);}
2506 else {Valid2 =
true;}
2507 if (verbose) cout <<
"CS_Delaunay validities i & j: " << DNN->
Valid(jet_i) <<
" " << Valid2 << endl;
2508 }
while ( !DNN->
Valid(jet_i) || !Valid2);
2509 if (! recombine_with_beam) {
2511 if (verbose) cout <<
"CS_Delaunay call _do_ij_recomb: " << jet_i <<
" " << jet_j <<
" " << SmallestDij << endl;
2512 _do_ij_recombination_step(jet_i, jet_j, SmallestDij, nn);
2515 points.push_back(newpoint);
2517 if (verbose) cout <<
"CS_Delaunay call _do_iB_recomb: " << jet_i <<
" " << SmallestDij << endl;
2518 _do_iB_recombination_step(jet_i, SmallestDij);
2520 if (
i == n-1) {
break;}
2521 vector<int> updated_neighbours;
2522 if (! recombine_with_beam) {
2525 points[points.size()-1], point3,
2526 updated_neighbours);
2527 if (static_cast<unsigned int> (point3) != points.size()-1) {
2528 throw Error(
"INTERNAL ERROR: point3 != points.size()-1");}
2532 vector<int>::iterator
it = updated_neighbours.begin();
2533 for (; it != updated_neighbours.end(); ++
it) {
2535 _add_ktdistance_to_map(ii, DijMap, DNN.
get());
2543 double yiB = jet_scale_for_algorithm(
_jets[ii]);
2548 if (DeltaR2 > 1.0) {
2551 double kt2i = jet_scale_for_algorithm(
_jets[ii]);
2553 if (kt2i <= jet_scale_for_algorithm(
_jets[jj])) {
2554 double dij = DeltaR2 * kt2i;
2566 using namespace std;
2568 vector<PseudoJet *> jetsp(
_jets.size());
2569 vector<int> indices(
_jets.size());
2570 for (
size_t i = 0;
i<
_jets.size();
i++) {
2574 for (
int n = jetsp.size();
n > 0;
n--) {
2576 double ymin = jet_scale_for_algorithm(*(jetsp[0]));
2578 for (
int i = 0;
i <
n;
i++) {
2579 double yiB = jet_scale_for_algorithm(*(jetsp[
i]));
2581 ymin = yiB; ii =
i; jj = -2;}
2583 for (
int i = 0;
i < n-1;
i++) {
2584 for (
int j =
i+1;
j <
n;
j++) {
2586 double y =
min(jet_scale_for_algorithm(*(jetsp[
i])),
2587 jet_scale_for_algorithm(*(jetsp[
j])))
2588 * jetsp[
i]->plain_distance(*jetsp[j])*
_invR2;
2589 if (y < ymin) {ymin =
y; ii =
i; jj =
j;}
2592 int newn = 2*jetsp.size() -
n;
2595 _do_ij_recombination_step(jetsp[ii]-&
_jets[0],
2596 jetsp[jj]-&
_jets[0], ymin, nn);
2597 jetsp[ii] = &
_jets[nn];
2598 jetsp[jj] = jetsp[n-1];
2600 indices[jj] = indices[n-1];
2602 _do_iB_recombination_step(jetsp[ii]-&
_jets[0], ymin);
2603 jetsp[ii] = jetsp[n-1];
2604 indices[ii] = indices[n-1];
2611 using namespace std;
2613 EEBriefJet *
const jetA,
const int _jets_index)
const {
2614 double E =
_jets[_jets_index].E();
2616 double p = jet_def().extra_param();
2617 switch (_jet_algorithm) {
2622 if (p <= 0 && scale < 1
e-300) scale = 1
e-300;
2623 scale = pow(scale,p);
2626 throw Error(
"Unrecognised jet algorithm");
2629 double norm =
_jets[_jets_index].modp2();
2631 norm = 1.0/sqrt(norm);
2632 jetA->
nx = norm *
_jets[_jets_index].px();
2633 jetA->
ny = norm *
_jets[_jets_index].py();
2634 jetA->
nz = norm *
_jets[_jets_index].pz();
2650 - jeta->
nz*jetb->
nz;
2655 _simple_N2_cluster<BriefJet>();
2658 _simple_N2_cluster<EEBriefJet>();
2663 using namespace std;
2665 if (_associated_cs != NULL
2666 && _associated_cs->will_delete_self_when_unused()) {
2667 _associated_cs->signal_imminent_self_deletion();
2668 delete _associated_cs;
2672 return (_associated_cs != NULL);
2675 return _associated_cs;
2678 if (!_associated_cs)
2679 throw Error(
"you requested information about the internal structure of a jet, but its associated ClusterSequence has gone out of scope.");
2680 return _associated_cs;
2683 return validated_cs()->has_partner(reference, partner);
2686 return validated_cs()->has_child(reference, child);
2689 return validated_cs()->has_parents(reference, parent1, parent2);
2692 if ((!has_associated_cluster_sequence()) || (!jet.has_associated_cluster_sequence()))
2693 throw Error(
"you requested information about the internal structure of a jet, but it is not associated with a ClusterSequence or its associated ClusterSequence has gone out of scope.");
2694 if (reference.associated_cluster_sequence() != jet.associated_cluster_sequence())
return false;
2695 return validated_cs()->object_in_jet(reference, jet);
2698 if (!has_associated_cluster_sequence())
2699 throw Error(
"you requested information about the internal structure of a jet, but it is not associated with a ClusterSequence or its associated ClusterSequence has gone out of scope.");
2703 return validated_cs()->constituents(reference);
2706 if (!has_associated_cluster_sequence())
2707 throw Error(
"you requested information about the internal structure of a jet, but it is not associated with a ClusterSequence or its associated ClusterSequence has gone out of scope.");
2711 return validated_cs()->exclusive_subjets(reference, dcut);
2714 return validated_cs()->n_exclusive_subjets(reference, dcut);
2717 return validated_cs()->exclusive_subjets_up_to(reference, nsub);
2720 return validated_cs()->exclusive_subdmerge(reference, nsub);
2723 return validated_cs()->exclusive_subdmerge_max(reference, nsub);
2726 PseudoJet dummy1, dummy2;
2727 return has_parents(reference, dummy1, dummy2);
2731 vector<PseudoJet> res;
2732 if (has_parents(reference, j1, j2)){
2744 using namespace std;
2752 if (jet->
next != NULL) {
2757 double default_size = max(0.1,
_Rparam);
2781 for (
int idphi = -1; idphi <=+1; idphi++) {
2792 for (
int idphi = -1; idphi <= +1; idphi++) {
2815 const int _jets_index) {
2816 _bj_set_jetinfo<>(jet, _jets_index);
2817 jet->tile_index =
_tile_index(jet->eta, jet->phi);
2819 jet->previous = NULL;
2820 jet->next = tile->
head;
2825 for (vector<Tile>::const_iterator
tile =
_tiles.begin();
2827 cout <<
"Tile " <<
tile -
_tiles.begin()<<
" = ";
2829 for (
TiledJet * jetI =
tile->head; jetI != NULL; jetI = jetI->next) {
2830 list.push_back(jetI-briefjets);
2832 sort(list.begin(),list.end());
2833 for (
unsigned int i = 0;
i < list.size();
i++) {cout <<
" "<<list[
i];}
2838 vector<int> & tile_union,
int & n_near_tiles)
const {
2839 for (
Tile *
const * near_tile =
_tiles[tile_index].begin_tiles;
2840 near_tile !=
_tiles[tile_index].end_tiles; near_tile++){
2841 tile_union[n_near_tiles] = *near_tile - &
_tiles[0];
2846 const int tile_index,
2847 vector<int> & tile_union,
int & n_near_tiles) {
2848 for (
Tile ** near_tile =
_tiles[tile_index].begin_tiles;
2849 near_tile !=
_tiles[tile_index].end_tiles; near_tile++){
2850 if (! (*near_tile)->tagged) {
2851 (*near_tile)->tagged =
true;
2852 tile_union[n_near_tiles] = *near_tile - &
_tiles[0];
2861 TiledJet * jetA = briefjets, * jetB;
2864 vector<int> tile_union(3*n_tile_neighbours);
2865 for (
int i = 0;
i<
n;
i++) {
2871 vector<Tile>::const_iterator
tile;
2872 for (tile =
_tiles.begin(); tile !=
_tiles.end(); tile++) {
2873 for (jetA = tile->head; jetA != NULL; jetA = jetA->
next) {
2874 for (jetB = tile->head; jetB != jetA; jetB = jetB->
next) {
2876 if (dist < jetA->NN_dist) {jetA->
NN_dist =
dist; jetA->
NN = jetB;}
2877 if (dist < jetB->NN_dist) {jetB->
NN_dist =
dist; jetB->NN = jetA;}
2880 for (
Tile ** RTile = tile->RH_tiles; RTile != tile->end_tiles; RTile++) {
2881 for (jetA = tile->head; jetA != NULL; jetA = jetA->
next) {
2882 for (jetB = (*RTile)->head; jetB != NULL; jetB = jetB->next) {
2884 if (dist < jetA->NN_dist) {jetA->
NN_dist =
dist; jetA->
NN = jetB;}
2885 if (dist < jetB->NN_dist) {jetB->
NN_dist =
dist; jetB->NN = jetA;}
2890 double * diJ =
new double[
n];
2892 for (
int i = 0;
i <
n;
i++) {
2896 int history_location = n-1;
2897 while (tail != head) {
2898 double diJ_min = diJ[0];
2899 int diJ_min_jet = 0;
2900 for (
int i = 1;
i <
n;
i++) {
2901 if (diJ[
i] < diJ_min) {diJ_min_jet =
i; diJ_min = diJ[
i];}
2904 jetA = & briefjets[diJ_min_jet];
2908 if (jetA < jetB) {
std::swap(jetA,jetB);}
2910 _do_ij_recombination_step(jetA->
_jets_index, jetB->_jets_index, diJ_min, nn);
2916 _do_iB_recombination_step(jetA->
_jets_index, diJ_min);
2919 int n_near_tiles = 0;
2922 bool sort_it =
false;
2934 sort(tile_union.begin(), tile_union.begin()+n_near_tiles);
2937 for (
int i = 1;
i < n_near_tiles;
i++) {
2938 if (tile_union[
i] != tile_union[nnn-1]) {
2939 tile_union[nnn] = tile_union[
i];
2950 diJ[jetA - head] = diJ[tail-head];
2958 for (
int itile = 0; itile < n_near_tiles; itile++) {
2959 Tile * tile_ptr = &
_tiles[tile_union[itile]];
2962 if (jetI->NN == jetA || (jetI->NN == jetB && jetB != NULL)) {
2963 jetI->NN_dist =
_R2;
2967 near_tile != tile_ptr->
end_tiles; near_tile++) {
2969 for (
TiledJet * jetJ = (*near_tile)->head;
2970 jetJ != NULL; jetJ = jetJ->
next) {
2972 if (dist < jetI->NN_dist && jetJ != jetI) {
2973 jetI->NN_dist =
dist; jetI->NN = jetJ;
2977 diJ[jetI-head] =
_bj_diJ(jetI);
2983 if (dist < jetI->NN_dist) {
2985 jetI->NN_dist =
dist;
2987 diJ[jetI-head] =
_bj_diJ(jetI);
2990 if (dist < jetB->NN_dist) {
2992 jetB->NN_dist =
dist;
2998 if (jetB != NULL) {diJ[jetB-head] =
_bj_diJ(jetB);}
3001 for (
TiledJet * jetJ = (*near_tile)->head;
3002 jetJ != NULL; jetJ = jetJ->
next) {
3003 if (jetJ->NN == tail) {jetJ->NN = jetA;}
3006 if (jetB != NULL) {diJ[jetB-head] =
_bj_diJ(jetB);}
3015 TiledJet * jetA = briefjets, * jetB;
3018 vector<int> tile_union(3*n_tile_neighbours);
3019 for (
int i = 0;
i<
n;
i++) {
3024 vector<Tile>::const_iterator
tile;
3025 for (tile =
_tiles.begin(); tile !=
_tiles.end(); tile++) {
3026 for (jetA = tile->head; jetA != NULL; jetA = jetA->
next) {
3027 for (jetB = tile->head; jetB != jetA; jetB = jetB->
next) {
3029 if (dist < jetA->NN_dist) {jetA->
NN_dist =
dist; jetA->
NN = jetB;}
3030 if (dist < jetB->NN_dist) {jetB->
NN_dist =
dist; jetB->NN = jetA;}
3033 for (
Tile ** RTile = tile->RH_tiles; RTile != tile->end_tiles; RTile++) {
3034 for (jetA = tile->head; jetA != NULL; jetA = jetA->
next) {
3035 for (jetB = (*RTile)->head; jetB != NULL; jetB = jetB->next) {
3037 if (dist < jetA->NN_dist) {jetA->
NN_dist =
dist; jetA->
NN = jetB;}
3038 if (dist < jetB->NN_dist) {jetB->
NN_dist =
dist; jetB->NN = jetA;}
3043 struct diJ_plus_link {
3047 diJ_plus_link * diJ =
new diJ_plus_link[
n];
3049 for (
int i = 0;
i <
n;
i++) {
3055 int history_location = n-1;
3057 diJ_plus_link * best, *stop;
3058 double diJ_min = diJ[0].diJ;
3061 for (diJ_plus_link * here = diJ+1; here != stop; here++) {
3062 if (here->diJ < diJ_min) {best = here; diJ_min = here->diJ;}
3069 if (jetA < jetB) {
std::swap(jetA,jetB);}
3071 _do_ij_recombination_step(jetA->
_jets_index, jetB->_jets_index, diJ_min, nn);
3077 _do_iB_recombination_step(jetA->
_jets_index, diJ_min);
3080 int n_near_tiles = 0;
3082 tile_union, n_near_tiles);
3086 tile_union,n_near_tiles);
3091 tile_union,n_near_tiles);
3097 for (
int itile = 0; itile < n_near_tiles; itile++) {
3098 Tile * tile_ptr = &
_tiles[tile_union[itile]];
3099 tile_ptr->
tagged =
false;
3102 if (jetI->NN == jetA || (jetI->NN == jetB && jetB != NULL)) {
3103 jetI->NN_dist =
_R2;
3107 near_tile != tile_ptr->
end_tiles; near_tile++) {
3109 for (
TiledJet * jetJ = (*near_tile)->head;
3110 jetJ != NULL; jetJ = jetJ->
next) {
3112 if (dist < jetI->NN_dist && jetJ != jetI) {
3113 jetI->NN_dist =
dist; jetI->NN = jetJ;
3117 diJ[jetI->diJ_posn].diJ =
_bj_diJ(jetI);
3124 if (dist < jetI->NN_dist) {
3126 jetI->NN_dist =
dist;
3128 diJ[jetI->diJ_posn].diJ =
_bj_diJ(jetI);
3131 if (dist < jetB->NN_dist) {
3133 jetB->NN_dist =
dist;
3139 if (jetB != NULL) {diJ[jetB->diJ_posn].diJ =
_bj_diJ(jetB);}
3148 TiledJet * jetA = briefjets, * jetB;
3151 vector<int> tile_union(3*n_tile_neighbours);
3152 for (
int i = 0;
i<
n;
i++) {
3157 vector<Tile>::const_iterator
tile;
3158 for (tile =
_tiles.begin(); tile !=
_tiles.end(); tile++) {
3159 for (jetA = tile->head; jetA != NULL; jetA = jetA->
next) {
3160 for (jetB = tile->head; jetB != jetA; jetB = jetB->
next) {
3162 if (dist < jetA->NN_dist) {jetA->
NN_dist =
dist; jetA->
NN = jetB;}
3163 if (dist < jetB->NN_dist) {jetB->
NN_dist =
dist; jetB->NN = jetA;}
3166 for (
Tile ** RTile = tile->RH_tiles; RTile != tile->end_tiles; RTile++) {
3167 for (jetA = tile->head; jetA != NULL; jetA = jetA->
next) {
3168 for (jetB = (*RTile)->head; jetB != NULL; jetB = jetB->next) {
3170 if (dist < jetA->NN_dist) {jetA->
NN_dist =
dist; jetA->
NN = jetB;}
3171 if (dist < jetB->NN_dist) {jetB->
NN_dist =
dist; jetB->NN = jetA;}
3176 vector<double> diJs(n);
3177 for (
int i = 0;
i <
n;
i++) {
3182 vector<TiledJet *> jets_for_minheap;
3183 jets_for_minheap.reserve(n);
3184 int history_location = n-1;
3187 jetA = head + minheap.
minloc();
3191 if (jetA < jetB) {
std::swap(jetA,jetB);}
3193 _do_ij_recombination_step(jetA->
_jets_index, jetB->_jets_index, diJ_min, nn);
3199 _do_iB_recombination_step(jetA->
_jets_index, diJ_min);
3202 minheap.
remove(jetA-head);
3203 int n_near_tiles = 0;
3205 tile_union, n_near_tiles);
3209 tile_union,n_near_tiles);
3221 tile_union,n_near_tiles);
3223 jetB->label_minheap_update_needed();
3224 jets_for_minheap.push_back(jetB);
3226 for (
int itile = 0; itile < n_near_tiles; itile++) {
3227 Tile * tile_ptr = &
_tiles[tile_union[itile]];
3228 tile_ptr->
tagged =
false;
3231 if (jetI->NN == jetA || (jetI->NN == jetB && jetB != NULL)) {
3232 jetI->NN_dist =
_R2;
3235 if (!jetI->minheap_update_needed()) {
3236 jetI->label_minheap_update_needed();
3237 jets_for_minheap.push_back(jetI);}
3240 near_tile != tile_ptr->
end_tiles; near_tile++) {
3242 for (
TiledJet * jetJ = (*near_tile)->head;
3243 jetJ != NULL; jetJ = jetJ->
next) {
3245 if (dist < jetI->NN_dist && jetJ != jetI) {
3246 jetI->NN_dist =
dist; jetI->NN = jetJ;
3256 if (dist < jetI->NN_dist) {
3258 jetI->NN_dist =
dist;
3261 if (!jetI->minheap_update_needed()) {
3262 jetI->label_minheap_update_needed();
3263 jets_for_minheap.push_back(jetI);}
3266 if (dist < jetB->NN_dist) {
3268 jetB->NN_dist =
dist;
3274 while (jets_for_minheap.size() > 0) {
3275 TiledJet * jetI = jets_for_minheap.back();
3276 jets_for_minheap.pop_back();
3286 using namespace std;
3289 : _pieces(initial_pieces){
3294 string str =
"Composite PseudoJet";
3301 vector<PseudoJet> all_constituents;
3302 for (
unsigned i = 0;
i <
_pieces.size();
i++) {
3304 vector<PseudoJet> constits =
_pieces[
i].constituents();
3305 copy(constits.begin(), constits.end(), back_inserter(all_constituents));
3307 all_constituents.push_back(
_pieces[
i]);
3310 return all_constituents;
3318 using namespace std;
3322 #if (!defined(FJCORE_HAVE_EXECINFO_H)) || defined(__FJCORE__)
3326 _message = message_in;
3327 if (_print_errors && _default_ostr){
3329 oss <<
"fjcore::Error: "<< message_in << endl;
3330 *_default_ostr << oss.str();
3331 _default_ostr->flush();
3335 #if (!defined(FJCORE_HAVE_EXECINFO_H)) || defined(__FJCORE__)
3337 _execinfo_undefined.warn(
"Error::set_print_backtrace(true) will not work with this build of FastJet");
3340 _print_backtrace = enabled;
3345 using namespace std;
3350 using namespace std;
3357 _jet_algorithm(jet_algorithm_in), _Rparam(R_in), _strategy(strategy_in) {
3361 if (R_in > max_allowable_R) {
3363 oss <<
"Requested R = " << R_in <<
" for jet definition is larger than max_allowable_R = " << max_allowable_R;
3364 throw Error(oss.str());
3367 unsigned int nparameters_expected = n_parameters_for_algorithm(jet_algorithm_in);
3368 if (nparameters != (
int) nparameters_expected){
3370 oss <<
"The jet algorithm you requested ("
3371 << jet_algorithm_in <<
") should be constructed with " << nparameters_expected
3372 <<
" parameter(s) but was called with " << nparameters <<
" parameter(s)\n";
3373 throw Error(oss.str());
3377 set_recombination_scheme(recomb_scheme_in);
3378 set_extra_param(0.0);
3382 return plugin()->is_spherical();
3391 name << description_no_recombiner();
3395 if (n_parameters_for_algorithm(jet_algorithm()) == 0)
3399 name << recombiner()->description();
3405 return plugin()->description();
3407 return "uninitialised JetDefinition (jet_algorithm=undefined_jet_algorithm)" ;
3409 name << algorithm_description(jet_algorithm());
3410 switch (n_parameters_for_algorithm(jet_algorithm())){
3411 case 0: name <<
" (NB: no R)";
break;
3412 case 1: name <<
" with R = " <<
R();
break;
3414 name <<
" with R = " <<
R();
3416 name <<
"and a special hack whereby particles with kt < "
3417 << extra_param() <<
"are treated as passive ghosts";
3419 name <<
", p = " << extra_param();
3428 case kt_algorithm:
return "Longitudinally invariant kt algorithm";
3430 case antikt_algorithm:
return "Longitudinally invariant anti-kt algorithm";
3431 case genkt_algorithm:
return "Longitudinally invariant generalised kt algorithm";
3437 throw Error(
"JetDefinition::algorithm_description(): unrecognized jet_algorithm");
3451 if (_shared_recombiner) _shared_recombiner.reset();
3455 assert(other_jet_def._recombiner ||
3457 if (other_jet_def._recombiner == 0){
3458 set_recombination_scheme(other_jet_def.recombination_scheme());
3461 _recombiner = other_jet_def._recombiner;
3463 _shared_recombiner.reset(other_jet_def._shared_recombiner);
3467 if (other_jd.recombination_scheme() != scheme)
return false;
3469 || (recombiner() == other_jd.recombiner());
3472 if (_recombiner == 0){
3473 throw Error(
"tried to call JetDefinition::delete_recombiner_when_unused() for a JetDefinition without a user-defined recombination scheme");
3474 }
else if (_shared_recombiner.get()) {
3475 throw Error(
"Error in JetDefinition::delete_recombiner_when_unused: the recombiner is already scheduled for deletion when unused (or was already set as shared)");
3477 _shared_recombiner.reset(_recombiner);
3481 throw Error(
"tried to call JetDefinition::delete_plugin_when_unused() for a JetDefinition without a plugin");
3483 _plugin_shared.reset(_plugin);
3488 return "E scheme recombination";
3490 return "pt scheme recombination";
3492 return "pt2 scheme recombination";
3494 return "Et scheme recombination";
3496 return "Et2 scheme recombination";
3498 return "boost-invariant pt scheme recombination";
3500 return "boost-invariant pt2 scheme recombination";
3502 return "pt-ordered Winner-Takes-All recombination";
3504 return "|3-momentum|-ordered Winner-Takes-All recombination";
3507 err <<
"DefaultRecombiner: unrecognized recombination scheme "
3509 throw Error(err.str());
3513 const PseudoJet & pa,
const PseudoJet & pb,
3514 PseudoJet & pab)
const {
3515 double weighta, weightb;
3516 switch(_recomb_scheme) {
3518 pab.reset(pa.px()+pb.px(),
3526 weighta = pa.perp();
3527 weightb = pb.perp();
3532 weighta = pa.perp2();
3533 weightb = pb.perp2();
3536 const PseudoJet & phard = (pa.pt2() >= pb.pt2()) ? pa : pb;
3537 pab.reset_PtYPhiM(pa.pt()+pb.pt(),
3538 phard.rap(), phard.phi(), phard.m());
3541 bool a_hardest = (pa.modp2() >= pb.modp2());
3542 const PseudoJet & phard = a_hardest ? pa : pb;
3543 const PseudoJet & psoft = a_hardest ? pb : pa;
3544 double modp_hard = phard.modp();
3545 double modp_ab = modp_hard + psoft.modp();
3546 if (phard.modp2()==0.0){
3547 pab.reset(0.0, 0.0, 0.0, phard.m());
3549 double scale = modp_ab/modp_hard;
3550 pab.reset(phard.px()*scale, phard.py()*scale, phard.pz()*scale,
3551 sqrt(modp_ab*modp_ab + phard.m2()));
3556 err <<
"DefaultRecombiner: unrecognized recombination scheme "
3558 throw Error(err.str());
3560 double perp_ab = pa.perp() + pb.perp();
3561 if (perp_ab != 0.0) {
3562 double y_ab = (weighta * pa.rap() + weightb * pb.rap())/(weighta+weightb);
3563 double phi_a = pa.phi(), phi_b = pb.phi();
3564 if (phi_a - phi_b >
pi) phi_b +=
twopi;
3565 if (phi_a - phi_b < -
pi) phi_b -=
twopi;
3566 double phi_ab = (weighta * phi_a + weightb * phi_b)/(weighta+weightb);
3567 pab.reset_PtYPhiM(perp_ab,y_ab,phi_ab);
3569 pab.reset(0.0, 0.0, 0.0, 0.0);
3573 switch(_recomb_scheme) {
3583 double newE = sqrt(p.perp2()+p.pz()*p.pz());
3584 p.reset_momentum(p.px(), p.py(), p.pz(), newE);
3590 double rescale = p.E()/sqrt(p.perp2()+p.pz()*p.pz());
3591 p.reset_momentum(rescale*p.px(), rescale*p.py(), rescale*p.pz(), p.E());
3596 err <<
"DefaultRecombiner: unrecognized recombination scheme "
3598 throw Error(err.str());
3602 throw Error(
"set_ghost_separation_scale not supported");
3606 if (pieces.size()>0){
3608 for (
unsigned int i=1;
i<pieces.size();
i++)
3615 PseudoJet
join(
const PseudoJet & j1,
3617 return join(vector<PseudoJet>(1,j1), recombiner);
3619 PseudoJet
join(
const PseudoJet & j1,
const PseudoJet & j2,
3621 vector<PseudoJet> pieces;
3622 pieces.push_back(j1);
3623 pieces.push_back(j2);
3624 return join(pieces, recombiner);
3626 PseudoJet
join(
const PseudoJet & j1,
const PseudoJet & j2,
const PseudoJet & j3,
3628 vector<PseudoJet> pieces;
3629 pieces.push_back(j1);
3630 pieces.push_back(j2);
3631 pieces.push_back(j3);
3632 return join(pieces, recombiner);
3634 PseudoJet
join(
const PseudoJet & j1,
const PseudoJet & j2,
const PseudoJet & j3,
const PseudoJet & j4,
3636 vector<PseudoJet> pieces;
3637 pieces.push_back(j1);
3638 pieces.push_back(j2);
3639 pieces.push_back(j3);
3640 pieces.push_back(j4);
3641 return join(pieces, recombiner);
3646 using namespace std;
3652 if (_this_warning_summary == 0) {
3653 _global_warnings_summary.push_back(
Summary(warning, 0));
3654 _this_warning_summary = & (_global_warnings_summary.back());
3656 if (_n_warn_so_far < _max_warn) {
3657 ostringstream warnstr;
3658 warnstr <<
"WARNING from FastJet: ";
3661 if (_n_warn_so_far == _max_warn) warnstr <<
" (LAST SUCH WARNING)";
3662 warnstr << std::endl;
3664 (*ostr) << warnstr.str();
3668 if (_this_warning_summary->second < numeric_limits<unsigned>::max()) {
3669 _this_warning_summary->second++;
3674 for (list<Summary>::const_iterator
it = _global_warnings_summary.begin();
3675 it != _global_warnings_summary.end();
it++) {
3676 str <<
it->second <<
" times: " <<
it->first << endl;
3685 using namespace std;
3687 for (
unsigned i = values.size();
i < _heap.size();
i++) {
3688 _heap[
i].value = std::numeric_limits<double>::max();
3689 _heap[
i].minloc = &(_heap[
i]);
3691 for (
unsigned i = 0;
i < values.size();
i++) {
3692 _heap[
i].value = values[
i];
3693 _heap[
i].minloc = &(_heap[
i]);
3695 for (
unsigned i = _heap.size()-1;
i > 0;
i--) {
3704 assert(loc < _heap.size());
3706 if (start->minloc != start && !(new_value < start->minloc->value)) {
3707 start->value = new_value;
3710 start->value = new_value;
3711 start->minloc =
start;
3712 bool change_made =
true;
3713 ValueLoc * heap_end = (&(_heap[0])) + _heap.size();
3714 while(change_made) {
3716 change_made =
false;
3717 if (here->
minloc == start) {
3718 here->
minloc = here; change_made =
true;
3720 ValueLoc * child = &(_heap[2*loc+1]);
3721 if (child < heap_end && child->minloc->value < here->
minloc->
value ) {
3723 change_made =
true;}
3725 if (child < heap_end && child->minloc->value < here->
minloc->
value ) {
3727 change_made =
true;}
3728 if (loc == 0) {
break;}
3740 using namespace std;
3741 PseudoJet::PseudoJet(
const double px_in,
const double py_in,
const double pz_in,
const double E_in) {
3746 this->_finish_init();
3750 _kt2 = this->px()*this->px() + this->py()*this->py();
3758 _phi = atan2(this->py(),this->px());
3760 if (_phi < 0.0) {_phi +=
twopi;}
3762 if (this->
E() == abs(this->pz()) && _kt2 == 0) {
3763 double MaxRapHere =
MaxRap + abs(this->pz());
3764 if (this->pz() >= 0.0) {_rap = MaxRapHere;}
else {_rap = -MaxRapHere;}
3766 double effective_m2 = max(0.0,
m2());
3767 double E_plus_pz = _E + abs(_pz);
3768 _rap = 0.5*log((_kt2 + effective_m2)/(E_plus_pz*E_plus_pz));
3769 if (_pz > 0) {_rap = - _rap;}
3773 valarray<double> mom(4);
3792 err <<
"PseudoJet subscripting: bad index (" << i <<
")";
3793 throw Error(err.str());
3798 if (px() == 0.0 && py() ==0.0)
return MaxRap;
3799 if (pz() == 0.0)
return 0.0;
3801 if (theta < 0) theta +=
pi;
3802 return -log(tan(theta/2));
3804 PseudoJet
operator+ (
const PseudoJet & jet1,
const PseudoJet & jet2) {
3805 return PseudoJet(jet1.px()+jet2.px(),
3806 jet1.py()+jet2.py(),
3807 jet1.pz()+jet2.pz(),
3808 jet1.E() +jet2.E() );
3810 PseudoJet
operator- (
const PseudoJet & jet1,
const PseudoJet & jet2) {
3811 return PseudoJet(jet1.px()-jet2.px(),
3812 jet1.py()-jet2.py(),
3813 jet1.pz()-jet2.pz(),
3814 jet1.E() -jet2.E() );
3817 jet._ensure_valid_rap_phi();
3818 PseudoJet coeff_times_jet(jet);
3819 coeff_times_jet *= coeff;
3820 return coeff_times_jet;
3826 return (1.0/coeff)*jet;
3829 _ensure_valid_rap_phi();
3837 (*this) *= 1.0/coeff;
3840 _px += other_jet._px;
3841 _py += other_jet._py;
3842 _pz += other_jet._pz;
3843 _E += other_jet._E ;
3847 _px -= other_jet._px;
3848 _py -= other_jet._py;
3849 _pz -= other_jet._pz;
3850 _E -= other_jet._E ;
3854 if (a.px() != b.px())
return false;
3855 if (a.py() != b.py())
return false;
3856 if (a.pz() != b.pz())
return false;
3857 if (a.E () != b.E ())
return false;
3858 if (a.user_index() != b.user_index())
return false;
3859 if (a.cluster_hist_index() != b.cluster_hist_index())
return false;
3860 if (a.user_info_ptr() != b.user_info_ptr())
return false;
3861 if (a.structure_ptr() != b.structure_ptr())
return false;
3866 throw Error(
"comparing a PseudoJet with a non-zero constant (double) is not allowed.");
3867 return (jet.px() == 0 && jet.py() == 0 &&
3868 jet.pz() == 0 && jet.E() == 0);
3871 if (prest.px() == 0.0 && prest.py() == 0.0 && prest.pz() == 0.0)
3873 double m_local = prest.m();
3875 double pf4 = ( px()*prest.px() + py()*prest.py()
3876 + pz()*prest.pz() +
E()*prest.E() )/m_local;
3877 double fn = (pf4 +
E()) / (prest.E() + m_local);
3878 _px += fn*prest.px();
3879 _py += fn*prest.py();
3880 _pz += fn*prest.pz();
3886 if (prest.px() == 0.0 && prest.py() == 0.0 && prest.pz() == 0.0)
3888 double m_local = prest.m();
3890 double pf4 = ( -px()*prest.px() - py()*prest.py()
3891 - pz()*prest.pz() +
E()*prest.E() )/m_local;
3892 double fn = (pf4 +
E()) / (prest.E() + m_local);
3893 _px -= fn*prest.px();
3894 _py -= fn*prest.py();
3895 _pz -= fn*prest.pz();
3901 return jeta.px() == jetb.px()
3902 && jeta.py() == jetb.py()
3903 && jeta.pz() == jetb.pz()
3904 && jeta.E() == jetb.E();
3907 _rap = rap_in; _phi =
phi_in;
3909 if (_phi < 0) _phi +=
twopi;
3913 double ptm = (m_in == 0) ? pt_in : sqrt(pt_in*pt_in+m_in*m_in);
3914 double exprap = exp(y_in);
3915 double pminus = ptm/exprap;
3916 double pplus = ptm*exprap;
3917 double px_local = pt_in*cos(phi_in);
3918 double py_local = pt_in*sin(phi_in);
3919 reset_momentum(px_local,py_local,0.5*(pplus-pminus),0.5*(pplus+pminus));
3920 set_cached_rap_phi(y_in,phi_in);
3924 double ptm = (m == 0) ? pt : sqrt(pt*pt+m*m);
3925 double exprap = exp(y);
3926 double pminus = ptm/exprap;
3927 double pplus = ptm*exprap;
3928 double px = pt*cos(phi);
3929 double py = pt*sin(phi);
3930 PseudoJet mom(px,py,0.5*(pplus-pminus),0.5*(pplus+pminus));
3931 mom.set_cached_rap_phi(y,phi);
3936 double dphi = abs(
phi() - other.phi());
3938 double drap = rap() - other.rap();
3939 distance = distance * (dphi*dphi + drap*drap);
3943 double dphi = abs(
phi() - other.phi());
3945 double drap = rap() - other.rap();
3946 return (dphi*dphi + drap*drap);
3949 double dphi = other.phi() -
phi();
3951 if (dphi < -
pi) dphi +=
twopi;
3956 return "standard PseudoJet (with no associated clustering information)";
3957 return _structure->description();
3960 return (_structure) && (_structure->has_associated_cluster_sequence());
3963 if (! has_associated_cluster_sequence())
return NULL;
3964 return _structure->associated_cluster_sequence();
3967 return (_structure) && (_structure->has_valid_cluster_sequence());
3970 return validated_structure_ptr()->validated_cs();
3973 _structure = structure_in;
3976 return bool(_structure);
3979 return _structure.get();
3982 return _structure.get();
3986 throw Error(
"Trying to access the structure of a PseudoJet which has no associated structure");
3987 return _structure.get();
3993 return validated_structure_ptr()->has_partner(*
this, partner);
3996 return validated_structure_ptr()->has_child(*
this, child);
3999 return validated_structure_ptr()->has_parents(*
this, parent1, parent2);
4002 return validated_structure_ptr()->object_in_jet(constituent, *
this);
4005 return validated_structure_ptr()->object_in_jet(*
this, jet);
4008 return (_structure) && (_structure->has_constituents());
4011 return validated_structure_ptr()->constituents(*
this);
4014 return (_structure) && (_structure->has_exclusive_subjets());
4017 return validated_structure_ptr()->exclusive_subjets(*
this, dcut);
4020 return validated_structure_ptr()->n_exclusive_subjets(*
this, dcut);
4023 return validated_structure_ptr()->exclusive_subjets_up_to(*
this, nsub);
4026 vector<PseudoJet> subjets = exclusive_subjets_up_to(nsub);
4027 if (
int(subjets.size()) < nsub) {
4029 err <<
"Requested " << nsub <<
" exclusive subjets, but there were only "
4030 << subjets.size() <<
" particles in the jet";
4031 throw Error(err.str());
4036 return validated_structure_ptr()->exclusive_subdmerge(*
this, nsub);
4039 return validated_structure_ptr()->exclusive_subdmerge_max(*
this, nsub);
4042 return ((_structure) && (_structure->has_pieces(*
this)));
4045 return validated_structure_ptr()->pieces(*
this);
4050 const vector<double> &
values) {
4052 sort(indices.begin(), indices.end(), index_sort_helper);
4055 vector<double> minus_kt2(jets.size());
4056 for (
size_t i = 0;
i < jets.size();
i++) {minus_kt2[
i] = -jets[
i].kt2();}
4060 vector<double> rapidities(jets.size());
4061 for (
size_t i = 0;
i < jets.size();
i++) {rapidities[
i] = jets[
i].rap();}
4065 vector<double> energies(jets.size());
4066 for (
size_t i = 0;
i < jets.size();
i++) {energies[
i] = -jets[
i].E();}
4070 vector<double>
pz(jets.size());
4071 for (
size_t i = 0;
i < jets.size();
i++) {
pz[
i] = jets[
i].pz();}
4076 for (
unsigned int i=0;
i<pieces.size();
i++)
4077 result += pieces[
i];
4083 return join(vector<PseudoJet>(1,j1));
4086 vector<PseudoJet>
pieces;
4088 pieces.push_back(j1);
4089 pieces.push_back(j2);
4090 return join(pieces);
4093 vector<PseudoJet>
pieces;
4095 pieces.push_back(j1);
4096 pieces.push_back(j2);
4097 pieces.push_back(j3);
4098 return join(pieces);
4101 vector<PseudoJet>
pieces;
4103 pieces.push_back(j1);
4104 pieces.push_back(j2);
4105 pieces.push_back(j3);
4106 pieces.push_back(j4);
4107 return join(pieces);
4110 using namespace std;
4116 throw Error(
"This PseudoJet structure is not associated with a valid ClusterSequence");
4119 throw Error(
"This PseudoJet structure has no implementation for has_partner");
4122 throw Error(
"This PseudoJet structure has no implementation for has_child");
4125 throw Error(
"This PseudoJet structure has no implementation for has_parents");
4128 throw Error(
"This PseudoJet structure has no implementation for is_inside");
4131 throw Error(
"This PseudoJet structure has no implementation for constituents");
4134 throw Error(
"This PseudoJet structure has no implementation for exclusive_subjets");
4137 throw Error(
"This PseudoJet structure has no implementation for n_exclusive_subjets");
4140 throw Error(
"This PseudoJet structure has no implementation for exclusive_subjets");
4143 throw Error(
"This PseudoJet structure has no implementation for exclusive_submerge");
4146 throw Error(
"This PseudoJet structure has no implementation for exclusive_submerge_max");
4149 throw Error(
"This PseudoJet structure has no implementation for pieces");
4153 #include <algorithm>
4154 using namespace std;
4157 std::vector<PseudoJet> result;
4158 const SelectorWorker * worker_local = validated_worker();
4159 if (worker_local->applies_jet_by_jet()) {
4160 for (std::vector<PseudoJet>::const_iterator jet = jets.begin();
4161 jet != jets.end(); jet++) {
4162 if (worker_local->pass(*jet)) result.push_back(*jet);
4165 std::vector<const PseudoJet *> jetptrs(jets.size());
4166 for (
unsigned i = 0;
i < jets.size();
i++) {
4167 jetptrs[
i] = & jets[
i];
4169 worker_local->terminator(jetptrs);
4170 for (
unsigned i = 0;
i < jetptrs.size();
i++) {
4171 if (jetptrs[
i]) result.push_back(jets[i]);
4178 const SelectorWorker * worker_local = validated_worker();
4179 if (worker_local->applies_jet_by_jet()) {
4180 for (
unsigned i = 0;
i < jets.size();
i++) {
4181 if (worker_local->pass(jets[
i])) n++;
4184 std::vector<const PseudoJet *> jetptrs(jets.size());
4185 for (
unsigned i = 0;
i < jets.size();
i++) {
4186 jetptrs[
i] = & jets[
i];
4188 worker_local->terminator(jetptrs);
4189 for (
unsigned i = 0;
i < jetptrs.size();
i++) {
4190 if (jetptrs[
i]) n++;
4197 const SelectorWorker * worker_local = validated_worker();
4198 if (worker_local->applies_jet_by_jet()) {
4199 for (
unsigned i = 0;
i < jets.size();
i++) {
4200 if (worker_local->pass(jets[
i])) this_sum += jets[i];
4203 std::vector<const PseudoJet *> jetptrs(jets.size());
4204 for (
unsigned i = 0;
i < jets.size();
i++) {
4205 jetptrs[
i] = & jets[
i];
4207 worker_local->terminator(jetptrs);
4208 for (
unsigned i = 0;
i < jetptrs.size();
i++) {
4209 if (jetptrs[
i]) this_sum += jets[
i];
4215 double this_sum = 0.0;
4216 const SelectorWorker * worker_local = validated_worker();
4217 if (worker_local->applies_jet_by_jet()) {
4218 for (
unsigned i = 0;
i < jets.size();
i++) {
4219 if (worker_local->pass(jets[
i])) this_sum += jets[i].
pt();
4222 std::vector<const PseudoJet *> jetptrs(jets.size());
4223 for (
unsigned i = 0;
i < jets.size();
i++) {
4224 jetptrs[
i] = & jets[
i];
4226 worker_local->terminator(jetptrs);
4227 for (
unsigned i = 0;
i < jetptrs.size();
i++) {
4228 if (jetptrs[
i]) this_sum += jets[
i].pt();
4234 std::vector<PseudoJet> & jets_that_pass,
4235 std::vector<PseudoJet> & jets_that_fail
4237 const SelectorWorker * worker_local = validated_worker();
4238 jets_that_pass.clear();
4239 jets_that_fail.clear();
4240 if (worker_local->applies_jet_by_jet()) {
4241 for (
unsigned i = 0;
i < jets.size();
i++) {
4242 if (worker_local->pass(jets[
i])) {
4243 jets_that_pass.push_back(jets[i]);
4245 jets_that_fail.push_back(jets[i]);
4249 std::vector<const PseudoJet *> jetptrs(jets.size());
4250 for (
unsigned i = 0;
i < jets.size();
i++) {
4251 jetptrs[
i] = & jets[
i];
4253 worker_local->terminator(jetptrs);
4254 for (
unsigned i = 0;
i < jetptrs.size();
i++) {
4256 jets_that_pass.push_back(jets[i]);
4258 jets_that_fail.push_back(jets[i]);
4264 if (! is_geometric())
return false;
4265 double rapmin, rapmax;
4266 get_rapidity_extent(rapmin, rapmax);
4267 return (rapmax != std::numeric_limits<double>::infinity())
4268 && (-rapmin != std::numeric_limits<double>::infinity());
4290 if (!applies_jet_by_jet())
4291 throw Error(
"Cannot apply this selector worker to an individual jet");
4292 return ! _s.pass(jet);
4295 virtual void terminator(vector<const PseudoJet *> & jets)
const {
4296 if (applies_jet_by_jet()){
4300 vector<const PseudoJet *> s_jets = jets;
4301 _s.worker()->terminator(s_jets);
4302 for (
unsigned int i=0;
i<s_jets.size();
i++){
4303 if (s_jets[
i]) jets[
i] = NULL;
4308 ostr <<
"!(" << _s.description() <<
")";
4323 _applies_jet_by_jet = _s1.applies_jet_by_jet() && _s2.applies_jet_by_jet();
4324 _takes_reference = _s1.takes_reference() || _s2.takes_reference();
4325 _is_geometric = _s1.is_geometric() && _s2.is_geometric();
4329 return _takes_reference;
4332 _s1.set_reference(centre);
4333 _s2.set_reference(centre);
4347 if (!applies_jet_by_jet())
4348 throw Error(
"Cannot apply this selector worker to an individual jet");
4349 return _s1.pass(jet) && _s2.pass(jet);
4351 virtual void terminator(vector<const PseudoJet *> & jets)
const {
4352 if (applies_jet_by_jet()){
4356 vector<const PseudoJet *> s1_jets = jets;
4357 _s1.worker()->terminator(s1_jets);
4358 _s2.worker()->terminator(jets);
4359 for (
unsigned int i=0;
i<jets.size();
i++){
4360 if (! s1_jets[
i]) jets[
i] = NULL;
4364 double s1min, s1max, s2min, s2max;
4365 _s1.get_rapidity_extent(s1min, s1max);
4366 _s2.get_rapidity_extent(s2min, s2max);
4367 rapmax =
min(s1max, s2max);
4368 rapmin = max(s1min, s2min);
4372 ostr <<
"(" << _s1.description() <<
" && " << _s2.description() <<
")";
4384 if (!applies_jet_by_jet())
4385 throw Error(
"Cannot apply this selector worker to an individual jet");
4386 return _s1.pass(jet) || _s2.pass(jet);
4389 return _s1.applies_jet_by_jet() && _s2.applies_jet_by_jet();
4391 virtual void terminator(vector<const PseudoJet *> & jets)
const {
4392 if (applies_jet_by_jet()){
4396 vector<const PseudoJet *> s1_jets = jets;
4397 _s1.worker()->terminator(s1_jets);
4398 _s2.worker()->terminator(jets);
4399 for (
unsigned int i=0;
i<jets.size();
i++){
4400 if (s1_jets[
i]) jets[
i] = s1_jets[
i];
4405 ostr <<
"(" << _s1.description() <<
" || " << _s2.description() <<
")";
4409 double s1min, s1max, s2min, s2max;
4410 _s1.get_rapidity_extent(s1min, s1max);
4411 _s2.get_rapidity_extent(s2min, s2max);
4412 rapmax = max(s1max, s2max);
4413 rapmin =
min(s1min, s2min);
4423 virtual void terminator(vector<const PseudoJet *> & jets)
const {
4424 if (applies_jet_by_jet()){
4428 _s2.worker()->terminator(jets);
4429 _s1.worker()->terminator(jets);
4433 ostr <<
"(" << _s1.description() <<
" * " << _s2.description() <<
")";
4459 template<
typename QuantityType>
4463 virtual bool pass(
const PseudoJet & jet)
const {
return _qmin(jet) >= _qmin.comparison_value();}
4466 ostr << _qmin.description() <<
" >= " << _qmin.description_value();
4473 template<
typename QuantityType>
4477 virtual bool pass(
const PseudoJet & jet)
const {
return _qmax(jet) <= _qmax.comparison_value();}
4480 ostr << _qmax.description() <<
" <= " << _qmax.description_value();
4487 template<
typename QuantityType>
4492 double q = _qmin(jet);
4493 return (q >= _qmin.comparison_value()) && (q <= _qmax.comparison_value());
4497 ostr << _qmin.description_value() <<
" <= " << _qmin.description() <<
" <= " << _qmax.description_value();
4576 rapmax = std::numeric_limits<double>::max();
4577 rapmin = _qmin.comparison_value();
4584 rapmax = _qmax.comparison_value();
4585 rapmin = -std::numeric_limits<double>::max();
4594 rapmax = _qmax.comparison_value();
4595 rapmin = _qmin.comparison_value();
4599 return twopi * (_qmax.comparison_value()-_qmin.comparison_value());
4622 rapmax = _qmax.comparison_value();
4623 rapmin = -_qmax.comparison_value();
4627 return twopi * 2 * _qmax.comparison_value();
4634 rapmax = _qmax.comparison_value();
4635 rapmin = -_qmax.comparison_value();
4639 return twopi * 2 * (_qmax.comparison_value()-max(_qmin.comparison_value(),0.0));
4688 _phispan = _phimax - _phimin;
4691 double dphi=jet.phi()-_phimin;
4693 if (dphi < 0) dphi +=
twopi;
4694 return (dphi <= _phispan);
4698 ostr << _phimin <<
" <= phi <= " << _phimax;
4714 _known_area = ((phimax-phimin >
twopi) ?
twopi : phimax-phimin) * (rapmax-rapmin);
4729 if (!applies_jet_by_jet())
4730 throw Error(
"Cannot apply this selector worker to an individual jet");
4733 virtual void terminator(vector<const PseudoJet *> & jets)
const {
4734 if (jets.size() < _n)
return;
4735 vector<double> minus_pt2(jets.size());
4736 vector<unsigned int> indices(jets.size());
4737 for (
unsigned int i=0;
i<jets.size();
i++){
4739 minus_pt2[
i] = jets[
i] ? -jets[
i]->perp2() : 0.0;
4742 partial_sort(indices.begin(), indices.begin()+_n, indices.end(), sort_helper);
4743 for (
unsigned int i=_n;
i<jets.size();
i++)
4744 jets[indices[
i]] = NULL;
4749 ostr << _n <<
" hardest";
4763 _is_initialised =
true;
4764 _reference = centre;
4775 if (! _is_initialised)
4776 throw Error(
"To use a SelectorCircle (or any selector that requires a reference), you first have to call set_reference(...)");
4777 return jet.squared_distance(_reference) <= _radius2;
4781 ostr <<
"distance from the centre <= " << sqrt(_radius2);
4785 if (! _is_initialised)
4786 throw Error(
"To use a SelectorCircle (or any selector that requires a reference), you first have to call set_reference(...)");
4787 rapmax = _reference.rap()+sqrt(_radius2);
4788 rapmin = _reference.rap()-sqrt(_radius2);
4794 return pi * _radius2;
4805 : _radius_in2(radius_in*radius_in), _radius_out2(radius_out*radius_out) {}
4808 if (! _is_initialised)
4809 throw Error(
"To use a SelectorDoughnut (or any selector that requires a reference), you first have to call set_reference(...)");
4810 double distance2 = jet.squared_distance(_reference);
4811 return (distance2 <= _radius_out2) && (distance2 >= _radius_in2);
4815 ostr << sqrt(_radius_in2) <<
" <= distance from the centre <= " << sqrt(_radius_out2);
4819 if (! _is_initialised)
4820 throw Error(
"To use a SelectorDoughnut (or any selector that requires a reference), you first have to call set_reference(...)");
4821 rapmax = _reference.rap()+sqrt(_radius_out2);
4822 rapmin = _reference.rap()-sqrt(_radius_out2);
4828 return pi * (_radius_out2-_radius_in2);
4841 if (! _is_initialised)
4842 throw Error(
"To use a SelectorStrip (or any selector that requires a reference), you first have to call set_reference(...)");
4843 return abs(jet.rap()-_reference.rap()) <= _delta;
4847 ostr <<
"|rap - rap_reference| <= " << _delta;
4851 if (! _is_initialised)
4852 throw Error(
"To use a SelectorStrip (or any selector that requires a reference), you first have to call set_reference(...)");
4853 rapmax = _reference.rap()+_delta;
4854 rapmin = _reference.rap()-_delta;
4860 return twopi * 2 * _delta;
4871 : _delta_rap(delta_rap), _delta_phi(delta_phi) {}
4874 if (! _is_initialised)
4875 throw Error(
"To use a SelectorRectangle (or any selector that requires a reference), you first have to call set_reference(...)");
4876 return (abs(jet.rap()-_reference.rap()) <= _delta_rap) && (abs(jet.delta_phi_to(_reference)) <= _delta_phi);
4880 ostr <<
"|rap - rap_reference| <= " << _delta_rap <<
" && |phi - phi_reference| <= " << _delta_phi ;
4884 if (! _is_initialised)
4885 throw Error(
"To use a SelectorRectangle (or any selector that requires a reference), you first have to call set_reference(...)");
4886 rapmax = _reference.rap()+_delta_rap;
4887 rapmin = _reference.rap()-_delta_rap;
4893 return 4 * _delta_rap * _delta_phi;
4906 if (! _is_initialised)
4907 throw Error(
"To use a SelectorPtFractionMin (or any selector that requires a reference), you first have to call set_reference(...)");
4908 return (jet.perp2() >= _fraction2*_reference.perp2());
4912 ostr <<
"pt >= " << sqrt(_fraction2) <<
"* pt_ref";
4933 _worker.reset(
new SW_And(*
this, b));
4937 _worker.reset(
new SW_Or(*
this, b));
4942 using namespace std;
4945 _cs(cs), _jets(cs.jets())
4950 #endif // INSTRUMENT2
4957 double default_size = max(0.1,
_Rparam)/2;
4961 #define _FJCORE_TILING25_USE_TILING_ANALYSIS_
4962 #ifdef _FASTJET_TILING25_USE_TILING_ANALYSIS_
4966 #else // not _FASTJET_TILING25_USE_TILING_ANALYSIS_
4969 const double maxrap = 7.0;
4970 for(
unsigned int i = 0;
i <
_jets.size();
i++) {
4971 double eta =
_jets[
i].rap();
4972 if (abs(eta) < maxrap) {
4977 #endif // _FASTJET_TILING25_USE_TILING_ANALYSIS_
4978 # define FJCORE_LAZY25_MIN3TILESY
4979 #ifdef FJCORE_LAZY25_MIN3TILESY
4986 #endif //FASTJET_LAZY25_MIN3TILESY
4991 #ifdef FJCORE_LAZY25_MIN3TILESY
4996 vector<bool> use_periodic_delta_phi(
_n_tiles_phi,
false);
4998 fill(use_periodic_delta_phi.begin(), use_periodic_delta_phi.end(),
true);
5000 use_periodic_delta_phi[0] =
true;
5001 use_periodic_delta_phi[1] =
true;
5010 tile->begin_tiles[0] =
tile;
5011 Tile25 ** pptile = & (tile->begin_tiles[0]);
5013 tile->surrounding_tiles = pptile;
5018 for (
int idphi = -2; idphi <=+2; idphi++) {
5027 for (
int idphi = -2; idphi <= +2; idphi++) {
5036 tile->RH_tiles = pptile;
5042 for (
int idphi = -2; idphi <= +2; idphi++) {
5048 for (
int idphi = -2; idphi <= +2; idphi++) {
5053 tile->end_tiles = pptile;
5054 tile->tagged =
false;
5055 tile->use_periodic_delta_phi = use_periodic_delta_phi[iphi];
5056 tile->max_NN_dist = 0;
5075 const int _jets_index) {
5076 _bj_set_jetinfo<>(jet, _jets_index);
5077 jet->tile_index =
_tile_index(jet->eta, jet->phi);
5079 jet->previous = NULL;
5080 jet->next = tile->head;
5081 if (jet->next != NULL) {jet->next->previous = jet;}
5087 tile->head = jet->
next;
5091 if (jet->
next != NULL) {
5096 for (vector<Tile25>::const_iterator
tile =
_tiles.begin();
5099 <<
" at " << setw(10) <<
tile->eta_centre <<
"," << setw(10) <<
tile->phi_centre
5102 for (
TiledJet * jetI =
tile->head; jetI != NULL; jetI = jetI->next) {
5103 list.push_back(jetI-briefjets);
5105 sort(list.begin(),list.end());
5106 for (
unsigned int i = 0;
i < list.size();
i++) {cout <<
" "<<list[
i];}
5111 vector<int> & tile_union,
int & n_near_tiles)
const {
5112 for (Tile25 *
const * near_tile =
_tiles[tile_index].begin_tiles;
5113 near_tile !=
_tiles[tile_index].end_tiles; near_tile++){
5114 tile_union[n_near_tiles] = *near_tile - &
_tiles[0];
5119 const int tile_index,
5120 vector<int> & tile_union,
int & n_near_tiles) {
5121 for (Tile25 ** near_tile =
_tiles[tile_index].begin_tiles;
5122 near_tile !=
_tiles[tile_index].end_tiles; near_tile++){
5123 if (! (*near_tile)->tagged) {
5124 (*near_tile)->tagged =
true;
5125 tile_union[n_near_tiles] = *near_tile - &
_tiles[0];
5132 vector<int> & tile_union,
int & n_near_tiles) {
5134 for (Tile25 ** near_tile = tile.begin_tiles; near_tile != tile.end_tiles; near_tile++){
5135 if ((*near_tile)->tagged)
continue;
5137 if (dist > (*near_tile)->max_NN_dist)
continue;
5138 (*near_tile)->tagged =
true;
5139 tile_union[n_near_tiles] = *near_tile - &
_tiles[0];
5149 #endif // INSTRUMENT2
5151 if (_tiles[bj->
tile_index].eta_centre == tile->eta_centre) deta = 0;
5152 else deta = std::abs(bj->
eta - tile->eta_centre) - _tile_half_size_eta;
5153 double dphi = std::abs(bj->
phi - tile->phi_centre);
5155 dphi -= _tile_half_size_phi;
5156 if (dphi < 0) dphi = 0;
5157 return dphi*dphi + deta*deta;
5161 if (dist < jetI->NN_dist) {
5167 jets_for_minheap.push_back(jetI);
5171 if (dist < jetX->NN_dist) {
5178 vector<TiledJet *> & jets_for_minheap) {
5183 jets_for_minheap.push_back(jetI);}
5185 for (Tile25 ** near_tile = tile_ptr->begin_tiles;
5186 near_tile != tile_ptr->end_tiles; near_tile++) {
5188 for (
TiledJet * jetJ = (*near_tile)->head;
5189 jetJ != NULL; jetJ = jetJ->
next) {
5191 if (dist < jetI->NN_dist && jetJ != jetI) {
5201 TiledJet * jetA = briefjets, * jetB;
5203 vector<int> tile_union(3*25);
5204 for (
int i = 0;
i<
n;
i++) {
5209 vector<Tile25>::iterator
tile;
5210 for (tile =
_tiles.begin(); tile !=
_tiles.end(); tile++) {
5211 for (jetA = tile->head; jetA != NULL; jetA = jetA->
next) {
5212 for (jetB = tile->head; jetB != jetA; jetB = jetB->
next) {
5214 if (dist < jetA->NN_dist) {jetA->
NN_dist =
dist; jetA->
NN = jetB;}
5215 if (dist < jetB->NN_dist) {jetB->
NN_dist =
dist; jetB->NN = jetA;}
5218 for (jetA = tile->head; jetA != NULL; jetA = jetA->
next) {
5219 if (jetA->
NN_dist > tile->max_NN_dist) tile->max_NN_dist = jetA->
NN_dist;
5222 for (tile =
_tiles.begin(); tile !=
_tiles.end(); tile++) {
5223 if (tile->use_periodic_delta_phi) {
5224 for (Tile25 ** RTile = tile->RH_tiles; RTile != tile->end_tiles; RTile++) {
5225 for (jetA = tile->head; jetA != NULL; jetA = jetA->
next) {
5227 bool relevant_for_jetA = dist_to_tile <= jetA->
NN_dist;
5228 bool relevant_for_RTile = dist_to_tile <= (*RTile)->max_NN_dist;
5229 if (relevant_for_jetA || relevant_for_RTile) {
5230 for (jetB = (*RTile)->head; jetB != NULL; jetB = jetB->next) {
5232 if (dist < jetA->NN_dist) {jetA->
NN_dist =
dist; jetA->
NN = jetB;}
5233 if (dist < jetB->NN_dist) {jetB->
NN_dist =
dist; jetB->NN = jetA;}
5239 for (Tile25 ** RTile = tile->RH_tiles; RTile != tile->end_tiles; RTile++) {
5240 for (jetA = tile->head; jetA != NULL; jetA = jetA->
next) {
5242 bool relevant_for_jetA = dist_to_tile <= jetA->
NN_dist;
5243 bool relevant_for_RTile = dist_to_tile <= (*RTile)->max_NN_dist;
5244 if (relevant_for_jetA || relevant_for_RTile) {
5245 for (jetB = (*RTile)->head; jetB != NULL; jetB = jetB->next) {
5247 if (dist < jetA->NN_dist) {jetA->
NN_dist =
dist; jetA->
NN = jetB;}
5248 if (dist < jetB->NN_dist) {jetB->
NN_dist =
dist; jetB->NN = jetA;}
5255 for (tile =
_tiles.begin(); tile !=
_tiles.end(); tile++) {
5256 tile->max_NN_dist = 0;
5257 for (jetA = tile->head; jetA != NULL; jetA = jetA->
next) {
5258 if (jetA->
NN_dist > tile->max_NN_dist) tile->max_NN_dist = jetA->
NN_dist;
5262 cout <<
"intermediate ncall, dtt = " << _ncall <<
" " << _ncall_dtt << endl;
5263 #endif // INSTRUMENT2
5264 vector<double> diJs(n);
5265 for (
int i = 0;
i <
n;
i++) {
5270 vector<TiledJet *> jets_for_minheap;
5271 jets_for_minheap.reserve(n);
5272 int history_location = n-1;
5275 jetA = head + minheap.
minloc();
5279 if (jetA < jetB) {
std::swap(jetA,jetB);}
5281 _cs.plugin_record_ij_recombination(jetA->
_jets_index, jetB->_jets_index, diJ_min, nn);
5290 minheap.
remove(jetA-head);
5291 int n_near_tiles = 0;
5293 Tile25 & jetB_tile =
_tiles[jetB->tile_index];
5294 for (Tile25 ** near_tile = jetB_tile.begin_tiles;
5295 near_tile != jetB_tile.end_tiles; near_tile++) {
5297 bool relevant_for_jetB = dist_to_tile <= jetB->NN_dist;
5298 bool relevant_for_near_tile = dist_to_tile <= (*near_tile)->max_NN_dist;
5299 bool relevant = relevant_for_jetB || relevant_for_near_tile;
5300 if (! relevant)
continue;
5301 tile_union[n_near_tiles] = *near_tile - &
_tiles[0];
5302 (*near_tile)->tagged =
true;
5304 for (
TiledJet * jetI = (*near_tile)->head; jetI != NULL; jetI = jetI->
next) {
5305 if (jetI->NN == jetA || jetI->
NN == jetB)
_set_NN(jetI, jets_for_minheap);
5310 int n_done_tiles = n_near_tiles;
5312 tile_union, n_near_tiles);
5315 tile_union,n_near_tiles);
5316 jetB->label_minheap_update_needed();
5317 jets_for_minheap.push_back(jetB);
5319 for (
int itile = 0; itile < n_done_tiles; itile++) {
5320 _tiles[tile_union[itile]].tagged =
false;
5322 for (
int itile = n_done_tiles; itile < n_near_tiles; itile++) {
5323 Tile25 * tile_ptr = &
_tiles[tile_union[itile]];
5324 tile_ptr->tagged =
false;
5325 for (
TiledJet * jetI = tile_ptr->head; jetI != NULL; jetI = jetI->
next) {
5326 if (jetI->NN == jetA || (jetI->NN == jetB && jetB != NULL)) {
5327 _set_NN(jetI, jets_for_minheap);
5331 while (jets_for_minheap.size() > 0) {
5332 TiledJet * jetI = jets_for_minheap.back();
5333 jets_for_minheap.pop_back();
5337 if (tile_I.max_NN_dist < jetI->
NN_dist) tile_I.max_NN_dist = jetI->
NN_dist;
5343 cout <<
"ncall, dtt = " << _ncall <<
" " << _ncall_dtt << endl;
5344 #endif // INSTRUMENT2
5350 using namespace std;
5351 #define _FJCORE_TILING2_USE_TILING_ANALYSIS_
5354 _cs(cs), _jets(cs.jets())
5359 #endif // INSTRUMENT2
5366 double default_size = max(0.1,
_Rparam);
5370 #ifdef _FJCORE_TILING2_USE_TILING_ANALYSIS_
5377 const double maxrap = 7.0;
5378 for(
unsigned int i = 0;
i <
_jets.size();
i++) {
5379 double eta =
_jets[
i].rap();
5380 if (abs(eta) < maxrap) {
5386 # define FJCORE_LAZY9_MIN2TILESY
5387 #ifdef FJCORE_LAZY9_MIN2TILESY
5394 #endif //FASTJET_LAZY9_MIN2TILESY
5399 #ifdef FJCORE_LAZY9_MIN2TILESY
5404 vector<bool> use_periodic_delta_phi(
_n_tiles_phi,
false);
5406 fill(use_periodic_delta_phi.begin(), use_periodic_delta_phi.end(),
true);
5408 use_periodic_delta_phi[0] =
true;
5424 for (
int idphi = -1; idphi <=+1; idphi++) {
5435 for (
int idphi = -1; idphi <= +1; idphi++) {
5462 const int _jets_index) {
5463 _bj_set_jetinfo<>(jet, _jets_index);
5464 jet->tile_index =
_tile_index(jet->eta, jet->phi);
5465 Tile2 * tile = &
_tiles[jet->tile_index];
5466 jet->previous = NULL;
5467 jet->next = tile->
head;
5478 if (jet->
next != NULL) {
5483 for (vector<Tile2>::const_iterator tile =
_tiles.begin();
5484 tile <
_tiles.end(); tile++) {
5485 cout <<
"Tile " << tile -
_tiles.begin()<<
" = ";
5487 for (
TiledJet * jetI = tile->head; jetI != NULL; jetI = jetI->
next) {
5488 list.push_back(jetI-briefjets);
5490 sort(list.begin(),list.end());
5491 for (
unsigned int i = 0;
i < list.size();
i++) {cout <<
" "<<list[
i];}
5496 vector<int> & tile_union,
int & n_near_tiles)
const {
5497 for (Tile2 *
const * near_tile =
_tiles[tile_index].begin_tiles;
5498 near_tile !=
_tiles[tile_index].end_tiles; near_tile++){
5499 tile_union[n_near_tiles] = *near_tile - &
_tiles[0];
5504 const int tile_index,
5505 vector<int> & tile_union,
int & n_near_tiles) {
5506 for (Tile2 ** near_tile =
_tiles[tile_index].begin_tiles;
5507 near_tile !=
_tiles[tile_index].end_tiles; near_tile++){
5508 if (! (*near_tile)->tagged) {
5509 (*near_tile)->tagged =
true;
5510 tile_union[n_near_tiles] = *near_tile - &
_tiles[0];
5517 vector<int> & tile_union,
int & n_near_tiles) {
5520 if ((*near_tile)->tagged)
continue;
5522 if (dist > (*near_tile)->max_NN_dist)
continue;
5523 (*near_tile)->tagged =
true;
5524 tile_union[n_near_tiles] = *near_tile - &
_tiles[0];
5534 #endif // INSTRUMENT2
5537 else deta = std::abs(bj->
eta - tile->
eta_centre) - _tile_half_size_eta;
5540 dphi -= _tile_half_size_phi;
5541 if (dphi < 0) dphi = 0;
5542 return dphi*dphi + deta*deta;
5546 if (dist < jetI->NN_dist) {
5552 jets_for_minheap.push_back(jetI);
5556 if (dist < jetX->NN_dist) {
5563 vector<TiledJet *> & jets_for_minheap) {
5568 jets_for_minheap.push_back(jetI);}
5571 near_tile != tile_ptr->
end_tiles; near_tile++) {
5573 for (
TiledJet * jetJ = (*near_tile)->head;
5574 jetJ != NULL; jetJ = jetJ->
next) {
5576 if (dist < jetI->NN_dist && jetJ != jetI) {
5586 TiledJet * jetA = briefjets, * jetB;
5588 vector<int> tile_union(3*n_tile_neighbours);
5589 for (
int i = 0;
i<
n;
i++) {
5594 vector<Tile2>::iterator
tile;
5595 for (tile =
_tiles.begin(); tile !=
_tiles.end(); tile++) {
5596 for (jetA = tile->head; jetA != NULL; jetA = jetA->
next) {
5597 for (jetB = tile->head; jetB != jetA; jetB = jetB->
next) {
5599 if (dist < jetA->NN_dist) {jetA->
NN_dist =
dist; jetA->
NN = jetB;}
5600 if (dist < jetB->NN_dist) {jetB->
NN_dist =
dist; jetB->NN = jetA;}
5603 for (jetA = tile->head; jetA != NULL; jetA = jetA->
next) {
5604 if (jetA->
NN_dist > tile->max_NN_dist) tile->max_NN_dist = jetA->
NN_dist;
5607 for (tile =
_tiles.begin(); tile !=
_tiles.end(); tile++) {
5608 if (tile->use_periodic_delta_phi) {
5609 for (Tile2 ** RTile = tile->RH_tiles; RTile != tile->end_tiles; RTile++) {
5610 for (jetA = tile->head; jetA != NULL; jetA = jetA->
next) {
5612 bool relevant_for_jetA = dist_to_tile <= jetA->
NN_dist;
5613 bool relevant_for_RTile = dist_to_tile <= (*RTile)->max_NN_dist;
5614 if (relevant_for_jetA || relevant_for_RTile) {
5615 for (jetB = (*RTile)->head; jetB != NULL; jetB = jetB->next) {
5617 if (dist < jetA->NN_dist) {jetA->
NN_dist =
dist; jetA->
NN = jetB;}
5618 if (dist < jetB->NN_dist) {jetB->
NN_dist =
dist; jetB->NN = jetA;}
5624 for (Tile2 ** RTile = tile->RH_tiles; RTile != tile->end_tiles; RTile++) {
5625 for (jetA = tile->head; jetA != NULL; jetA = jetA->
next) {
5627 bool relevant_for_jetA = dist_to_tile <= jetA->
NN_dist;
5628 bool relevant_for_RTile = dist_to_tile <= (*RTile)->max_NN_dist;
5629 if (relevant_for_jetA || relevant_for_RTile) {
5630 for (jetB = (*RTile)->head; jetB != NULL; jetB = jetB->next) {
5632 if (dist < jetA->NN_dist) {jetA->
NN_dist =
dist; jetA->
NN = jetB;}
5633 if (dist < jetB->NN_dist) {jetB->
NN_dist =
dist; jetB->NN = jetA;}
5640 for (tile =
_tiles.begin(); tile !=
_tiles.end(); tile++) {
5641 tile->max_NN_dist = 0;
5642 for (jetA = tile->head; jetA != NULL; jetA = jetA->
next) {
5643 if (jetA->
NN_dist > tile->max_NN_dist) tile->max_NN_dist = jetA->
NN_dist;
5647 cout <<
"intermediate ncall, dtt = " << _ncall <<
" " << _ncall_dtt << endl;
5648 #endif // INSTRUMENT2
5649 vector<double> diJs(n);
5650 for (
int i = 0;
i <
n;
i++) {
5655 vector<TiledJet *> jets_for_minheap;
5656 jets_for_minheap.reserve(n);
5657 int history_location = n-1;
5660 jetA = head + minheap.
minloc();
5664 if (jetA < jetB) {
std::swap(jetA,jetB);}
5666 _cs.plugin_record_ij_recombination(jetA->
_jets_index, jetB->_jets_index, diJ_min, nn);
5675 minheap.
remove(jetA-head);
5676 int n_near_tiles = 0;
5678 Tile2 & jetB_tile =
_tiles[jetB->tile_index];
5680 near_tile != jetB_tile.
end_tiles; near_tile++) {
5682 bool relevant_for_jetB = dist_to_tile <= jetB->NN_dist;
5683 bool relevant_for_near_tile = dist_to_tile <= (*near_tile)->max_NN_dist;
5684 bool relevant = relevant_for_jetB || relevant_for_near_tile;
5685 if (! relevant)
continue;
5686 tile_union[n_near_tiles] = *near_tile - &
_tiles[0];
5687 (*near_tile)->tagged =
true;
5689 for (
TiledJet * jetI = (*near_tile)->head; jetI != NULL; jetI = jetI->
next) {
5690 if (jetI->NN == jetA || jetI->
NN == jetB)
_set_NN(jetI, jets_for_minheap);
5695 int n_done_tiles = n_near_tiles;
5697 tile_union, n_near_tiles);
5700 tile_union,n_near_tiles);
5701 jetB->label_minheap_update_needed();
5702 jets_for_minheap.push_back(jetB);
5704 for (
int itile = 0; itile < n_done_tiles; itile++) {
5705 _tiles[tile_union[itile]].tagged =
false;
5707 for (
int itile = n_done_tiles; itile < n_near_tiles; itile++) {
5708 Tile2 * tile_ptr = &
_tiles[tile_union[itile]];
5709 tile_ptr->
tagged =
false;
5711 if (jetI->NN == jetA || (jetI->NN == jetB && jetB != NULL)) {
5712 _set_NN(jetI, jets_for_minheap);
5716 while (jets_for_minheap.size() > 0) {
5717 TiledJet * jetI = jets_for_minheap.back();
5718 jets_for_minheap.pop_back();
5728 cout <<
"ncall, dtt = " << _ncall <<
" " << _ncall_dtt << endl;
5729 #endif // INSTRUMENT2
5733 using namespace std;
5736 _cs(cs), _jets(cs.jets())
5744 double default_size = max(0.1,
_Rparam);
5750 const double maxrap = 7.0;
5751 for(
unsigned int i = 0;
i <
_jets.size();
i++) {
5752 double eta =
_jets[
i].rap();
5753 if (abs(eta) < maxrap) {
5764 vector<bool> use_periodic_delta_phi(
_n_tiles_phi,
false);
5766 fill(use_periodic_delta_phi.begin(), use_periodic_delta_phi.end(),
true);
5768 use_periodic_delta_phi[0] =
true;
5822 tile->
eta_max = (ieta+1)*_tile_size_eta;
5824 tile->
phi_max = (iphi+1)*_tile_size_phi;
5841 const int _jets_index) {
5842 _bj_set_jetinfo<>(jet, _jets_index);
5843 jet->tile_index =
_tile_index(jet->eta, jet->phi);
5845 jet->previous = NULL;
5846 jet->next = tile->
head;
5857 if (jet->
next != NULL) {
5862 for (vector<Tile>::const_iterator tile =
_tiles.begin();
5863 tile <
_tiles.end(); tile++) {
5864 cout <<
"Tile " << tile -
_tiles.begin()<<
" = ";
5867 list.push_back(jetI-briefjets);
5869 sort(list.begin(),list.end());
5870 for (
unsigned int i = 0;
i < list.size();
i++) {cout <<
" "<<list[
i];}
5875 vector<int> & tile_union,
int & n_near_tiles)
const {
5877 near_tile !=
_tiles[tile_index].end_tiles; near_tile++){
5878 tile_union[n_near_tiles] = near_tile->first - &
_tiles[0];
5883 const int tile_index,
5884 vector<int> & tile_union,
int & n_near_tiles) {
5886 near_tile !=
_tiles[tile_index].end_tiles; near_tile++){
5887 if (! (near_tile->first)->tagged) {
5888 (near_tile->first)->tagged =
true;
5889 tile_union[n_near_tiles] = near_tile->first - &
_tiles[0];
5896 vector<int> & tile_union,
int & n_near_tiles) {
5899 if ((near_tile->first)->tagged)
continue;
5900 double dist = (tile.*(near_tile->second))(jet) - tile_edge_security_margin;
5901 if (dist > (near_tile->first)->max_NN_dist)
continue;
5902 (near_tile->first)->tagged =
true;
5903 tile_union[n_near_tiles] = near_tile->first - &
_tiles[0];
5908 ostr <<
"j" << setw(3) << jet.
_jets_index <<
":pt2,rap,phi=" ; ostr.flush();
5909 ostr << jet.
kt2 <<
","; ostr.flush();
5910 ostr << jet.
eta <<
","; ostr.flush();
5911 ostr << jet.
phi; ostr.flush();
5912 ostr <<
", tile=" << jet.
tile_index; ostr.flush();
5917 if (dist < jetI->NN_dist) {
5923 jets_for_minheap.push_back(jetI);
5927 if (dist < jetX->NN_dist) {
5934 vector<TiledJet *> & jets_for_minheap) {
5939 jets_for_minheap.push_back(jetI);}
5942 near_tile != tile_ptr->
end_tiles; near_tile++) {
5943 if (jetI->
NN_dist < (tile_ptr->*(near_tile->second))(jetI))
continue;
5944 for (
TiledJet * jetJ = (near_tile->first)->head;
5945 jetJ != NULL; jetJ = jetJ->
next) {
5947 if (dist < jetI->NN_dist && jetJ != jetI) {
5956 TiledJet * jetA = briefjets, * jetB;
5958 vector<int> tile_union(3*n_tile_neighbours);
5959 for (
int i = 0;
i<
n;
i++) {
5964 vector<Tile>::iterator
tile;
5965 for (tile =
_tiles.begin(); tile !=
_tiles.end(); tile++) {
5966 for (jetA = tile->head; jetA != NULL; jetA = jetA->
next) {
5967 for (jetB = tile->head; jetB != jetA; jetB = jetB->
next) {
5969 if (dist < jetA->NN_dist) {jetA->
NN_dist =
dist; jetA->
NN = jetB;}
5970 if (dist < jetB->NN_dist) {jetB->
NN_dist =
dist; jetB->NN = jetA;}
5973 for (jetA = tile->head; jetA != NULL; jetA = jetA->
next) {
5974 if (jetA->
NN_dist > tile->max_NN_dist) tile->max_NN_dist = jetA->
NN_dist;
5977 for (tile =
_tiles.begin(); tile !=
_tiles.end(); tile++) {
5978 if (tile->use_periodic_delta_phi) {
5980 RTileFnPair != tile->end_tiles; RTileFnPair++) {
5981 Tile *RTile = RTileFnPair->first;
5982 for (jetA = tile->head; jetA != NULL; jetA = jetA->
next) {
5983 double dist_to_tile = ((*tile).*(RTileFnPair->second))(jetA);
5984 bool relevant_for_jetA = dist_to_tile <= jetA->
NN_dist;
5985 bool relevant_for_RTile = dist_to_tile <= RTile->
max_NN_dist;
5986 if (relevant_for_jetA || relevant_for_RTile) {
5987 for (jetB = RTile->
head; jetB != NULL; jetB = jetB->
next) {
5989 if (dist < jetA->NN_dist) {jetA->
NN_dist =
dist; jetA->
NN = jetB;}
5990 if (dist < jetB->NN_dist) {jetB->
NN_dist =
dist; jetB->NN = jetA;}
5997 RTileFnPair != tile->end_tiles; RTileFnPair++) {
5998 Tile *RTile = RTileFnPair->first;
5999 for (jetA = tile->head; jetA != NULL; jetA = jetA->
next) {
6000 double dist_to_tile = ((*tile).*(RTileFnPair->second))(jetA);
6001 bool relevant_for_jetA = dist_to_tile <= jetA->
NN_dist;
6002 bool relevant_for_RTile = dist_to_tile <= RTile->
max_NN_dist;
6003 if (relevant_for_jetA || relevant_for_RTile) {
6004 for (jetB = RTile->
head; jetB != NULL; jetB = jetB->
next) {
6006 if (dist < jetA->NN_dist) {jetA->
NN_dist =
dist; jetA->
NN = jetB;}
6007 if (dist < jetB->NN_dist) {jetB->
NN_dist =
dist; jetB->NN = jetA;}
6014 for (tile =
_tiles.begin(); tile !=
_tiles.end(); tile++) {
6015 tile->max_NN_dist = 0;
6016 for (jetA = tile->head; jetA != NULL; jetA = jetA->
next) {
6017 if (jetA->
NN_dist > tile->max_NN_dist) tile->max_NN_dist = jetA->
NN_dist;
6020 vector<double> diJs(n);
6021 for (
int i = 0;
i <
n;
i++) {
6026 vector<TiledJet *> jets_for_minheap;
6027 jets_for_minheap.reserve(n);
6028 int history_location = n-1;
6031 jetA = head + minheap.
minloc();
6035 if (jetA < jetB) {
std::swap(jetA,jetB);}
6037 _cs.plugin_record_ij_recombination(jetA->
_jets_index, jetB->_jets_index, diJ_min, nn);
6046 minheap.
remove(jetA-head);
6047 int n_near_tiles = 0;
6049 tile_union, n_near_tiles);
6052 tile_union,n_near_tiles);
6053 jetB->label_minheap_update_needed();
6054 jets_for_minheap.push_back(jetB);
6059 near_tile_fn_pair != jetB_tile.
end_tiles; near_tile_fn_pair++) {
6060 Tile * near_tile = near_tile_fn_pair->first;
6061 double dist_to_tile = (jetB_tile.*(near_tile_fn_pair->second))(jetB);
6062 bool relevant_for_jetB = dist_to_tile <= jetB->NN_dist;
6063 bool relevant_for_near_tile = dist_to_tile <= near_tile->
max_NN_dist;
6064 bool relevant = relevant_for_jetB || relevant_for_near_tile;
6068 if (jetI->NN == jetA || jetI->
NN == jetB)
_set_NN(jetI, jets_for_minheap);
6071 near_tile->
tagged =
false;
6081 for (
int itile = 0; itile < n_near_tiles; itile++) {
6082 Tile * tile_ptr = &
_tiles[tile_union[itile]];
6083 if (!tile_ptr->
tagged)
continue;
6084 tile_ptr->
tagged =
false;
6086 if (jetI->NN == jetA || (jetI->NN == jetB && jetB != NULL)) {
6087 _set_NN(jetI, jets_for_minheap);
6091 while (jets_for_minheap.size() > 0) {
6092 TiledJet * jetI = jets_for_minheap.back();
6093 jets_for_minheap.pop_back();
6107 using namespace std;
6110 _determine_rapidity_extent(cs.jets());
6113 _determine_rapidity_extent(particles);
6118 vector<double> counts(nbins, 0);
6119 _minrap = numeric_limits<double>::max();
6120 _maxrap = -numeric_limits<double>::max();
6122 for (
unsigned i = 0;
i < particles.size();
i++) {
6123 if (particles[
i].
E() == abs(particles[
i].pz()))
continue;
6124 double rap = particles[
i].rap();
6125 if (rap < _minrap) _minrap = rap;
6126 if (rap > _maxrap) _maxrap = rap;
6127 ibin = int(rap+nrap);
6128 if (ibin < 0) ibin = 0;
6129 if (ibin >= nbins) ibin = nbins - 1;
6132 double max_in_bin = 0;
6133 for (ibin = 0; ibin <
nbins; ibin++) {
6134 if (max_in_bin < counts[ibin]) max_in_bin = counts[
ibin];
6136 const double allowed_max_fraction = 0.25;
6137 const double min_multiplicity = 4;
6138 double allowed_max_cumul = floor(max(max_in_bin * allowed_max_fraction, min_multiplicity));
6139 if (allowed_max_cumul > max_in_bin) allowed_max_cumul = max_in_bin;
6140 double cumul_lo = 0;
6142 for (ibin = 0; ibin <
nbins; ibin++) {
6143 cumul_lo += counts[
ibin];
6144 if (cumul_lo >= allowed_max_cumul) {
6145 double y = ibin-nrap;
6146 if (y > _minrap) _minrap =
y;
6151 _cumul2 += cumul_lo*cumul_lo;
6153 double cumul_hi = 0;
6154 for (ibin = nbins-1; ibin >= 0; ibin--) {
6155 cumul_hi += counts[
ibin];
6156 if (cumul_hi >= allowed_max_cumul) {
6157 double y = ibin-nrap+1;
6158 if (y < _maxrap) _maxrap =
y;
6164 assert(ibin_hi >= ibin_lo);
6165 if (ibin_hi == ibin_lo) {
6166 _cumul2 = pow(
double(cumul_lo + cumul_hi - counts[ibin_hi]), 2);
6168 _cumul2 += cumul_hi*cumul_hi;
6169 for (ibin = ibin_lo+1; ibin < ibin_hi; ibin++) {
6170 _cumul2 += counts[
ibin]*counts[
ibin];