59 #include "protocones.h"
60 #include "siscone_error.h"
64 #include "circulator.h"
106 if (
hc!=NULL)
delete hc;
122 multiple_centre_done.clear();
156 for (p_idx=0;p_idx<
n_part;p_idx++){
177 }
while (!update_cone());
180 return proceed_with_stability();
198 int Cstable_cones::init_cone(){
206 prepare_cocircular_lists();
216 centre_idx = first_cone;
221 compute_cone_contents();
233 int Cstable_cones::test_cone(){
234 Creference weighted_cone_ref;
251 cone_candidate = cone;
256 cone_candidate = cone;
257 cone_candidate+= *
parent + *child;
261 cone_candidate = cone + *
parent;
262 hc->
insert(&cone_candidate, parent, child,
true,
false);
264 cone_candidate = cone + *child;
265 hc->
insert(&cone_candidate, parent, child,
false,
true);
279 int Cstable_cones::update_cone(){
284 if (centre_idx==first_cone)
299 dpt += fabs(child->px)+fabs(child->py);
311 if (cocircular_check())
312 return update_cone();
327 dpt += fabs(child->px)+fabs(child->py);
336 recompute_cone_contents();
352 int Cstable_cones::proceed_with_stability(){
357 for (i=0;i<=
hc->
mask;i++){
366 #ifdef USE_QUADTREE_FOR_STABILITY_TEST
368 if (quadtree->circle_intersect(elm->eta, elm->phi,
R2)==elm->ref){
371 if (circle_intersect(elm->eta, elm->phi)==elm->
ref){
378 protocones.push_back(Cmomentum(elm->eta, elm->phi, elm->ref));
390 #ifdef DEBUG_STABLE_CONES
392 nb_hash_occupied =
hc->n_occupied_cells;
412 bool cocircular_pt_less(Cmomentum *v1, Cmomentum *v2){
413 return v1->perp2() < v2->perp2();
424 void Cstable_cones::prepare_cocircular_lists() {
425 circulator<vector<Cvicinity_elm*>::iterator > here(
vicinity.begin(),
429 circulator<vector<Cvicinity_elm*>::iterator > search(here);
432 Cvicinity_elm* here_pntr = *here();
433 search.set_position(here);
439 if ( abs_dphi((*search())->angle, here_pntr->angle) <
440 here_pntr->cocircular_range
441 && search() != here()) {
442 (*search())->cocircular.push_back(here_pntr);
449 search.set_position(here);
452 if ( abs_dphi((*search())->angle, here_pntr->angle) <
453 here_pntr->cocircular_range
454 && search() != here()) {
455 (*search())->cocircular.push_back(here_pntr);
462 }
while (here() !=
vicinity.begin());
473 void Cstable_cones::test_cone_cocircular(Cmomentum & borderless_cone,
474 list<Cmomentum *> & border_list) {
475 vector<Cborder_store> border_vect;
477 border_vect.reserve(border_list.size());
478 for (list<Cmomentum *>::iterator it = border_list.begin();
479 it != border_list.end(); it++) {
480 border_vect.push_back(Cborder_store(*it, centre->
eta, centre->
phi));
484 sort(border_vect.begin(), border_vect.end());
488 circulator<vector<Cborder_store>::iterator >
489 start(border_vect.begin(), border_vect.begin(),border_vect.end());
490 circulator<vector<Cborder_store>::iterator > mid(start), end(start);
493 Cmomentum candidate = borderless_cone;
494 candidate.build_etaphi();
495 if (candidate.ref.not_empty())
496 test_stability(candidate, border_vect);
502 mid()->is_in =
false;
503 }
while (++mid != start);
506 candidate = borderless_cone;
507 while (++mid != start) {
510 candidate += *(mid()->mom);
511 test_stability(candidate, border_vect);
514 }
while (++start != end);
519 candidate += *(mid()->mom);
520 test_stability(candidate, border_vect);
530 void Cstable_cones::test_stability(Cmomentum & candidate,
const vector<Cborder_store> & border_vect) {
533 candidate.build_etaphi();
536 for (
unsigned i = 0; i < border_vect.size(); i++) {
537 if (is_inside(&candidate, border_vect[i].mom) ^ (border_vect[i].is_in)) {
553 bool Cstable_cones::cocircular_check(){
572 dpt += fabs(child->px)+fabs(child->py);
579 list<Cvicinity_inclusion *> removed_from_cone;
580 list<Cvicinity_inclusion *> put_in_border;
581 list<Cmomentum *> border_list;
583 Cmomentum cone_removal;
584 Cmomentum border = *
parent;
585 border_list.push_back(parent);
592 for(list<Cvicinity_elm *>::iterator it = centre->
cocircular.begin();
595 if ((*it)->is_inside->cone) {
596 cone_removal += *((*it)->v);
597 (*it)->is_inside->cone =
false;
598 removed_from_cone.push_back((*it)->is_inside);
605 if (!(*it)->is_inside->cocirc) {
606 border += *((*it)->v);
607 (*it)->is_inside->cocirc =
true;
608 put_in_border.push_back((*it)->is_inside);
609 border_list.push_back((*it)->v);
615 Cmomentum borderless_cone = cone;
616 borderless_cone -= cone_removal;
617 bool consider =
true;
618 for (
unsigned int i=0;i<multiple_centre_done.size();i++){
619 if ((multiple_centre_done[i].first ==borderless_cone.ref) &&
620 (multiple_centre_done[i].second==border.ref))
627 multiple_centre_done.push_back(pair<Creference,Creference>(borderless_cone.ref,
631 double local_dpt = fabs(cone_removal.px) + fabs(cone_removal.py);
632 double total_dpt = dpt + local_dpt;
634 recompute_cone_contents_if_needed(borderless_cone, total_dpt);
635 if (total_dpt == 0) {
638 cone = borderless_cone + cone_removal;
642 test_cone_cocircular(borderless_cone, border_list);
647 for(list<Cvicinity_inclusion *>::iterator is_in = removed_from_cone.begin();
648 is_in != removed_from_cone.end(); is_in++) {
649 (*is_in)->cone =
true;
653 for(list<Cvicinity_inclusion *>::iterator is_in = put_in_border.begin();
654 is_in != put_in_border.end(); is_in++) {
655 (*is_in)->cocirc =
false;
679 void Cstable_cones::compute_cone_contents() {
680 circulator<vector<Cvicinity_elm*>::iterator >
683 circulator<vector<Cvicinity_elm*>::iterator > here(start);
695 if (!(*here())->side) ((*here())->is_inside->cone) = 1;
701 if ((*here())->side) ((*here())->is_inside->cone) = 0;
702 }
while (here != start);
707 recompute_cone_contents();
718 void Cstable_cones::recompute_cone_contents(){
746 void Cstable_cones::recompute_cone_contents_if_needed(Cmomentum & this_cone,
749 if (this_dpt >
PT_TSHOLD*(fabs(this_cone.px)+fabs(this_cone.py))) {
751 this_cone = Cmomentum();
754 this_cone = Cmomentum();
791 Creference Cstable_cones::circle_intersect(
double cx,
double cy){
792 Creference intersection;
798 dx =
plist[i].eta - cx;
799 dy = fabs(
plist[i].phi - cy);
807 intersection+=
plist[i].ref;
821 inline bool Cstable_cones::is_inside(Cmomentum *centre_in, Cmomentum *v){
824 dx = centre_in->eta - v->eta;
825 dy = fabs(centre_in->phi - v->phi);
829 return dx*dx+dy*dy<
R2;
839 inline double abs_dangle(
double &angle1,
double &angle2){
842 dphi = fabs(angle1-angle2);