59 #include "protocones.h"
60 #include "siscone_error.h"
64 #include "circulator.h"
97 : Cvicinity(_particle_list){
105 Cstable_cones::~Cstable_cones(){
106 if (hc!=NULL)
delete hc;
114 void Cstable_cones::init(vector<Cmomentum> &_particle_list){
119 if (protocones.size()!=0)
122 multiple_centre_done.clear();
125 set_particle_list(_particle_list);
141 int Cstable_cones::get_stable_cones(
double _radius){
156 for (p_idx=0;p_idx<n_part;p_idx++){
159 build(&plist[p_idx], 2.0*R);
164 if (vicinity_size==0){
165 protocones.push_back(*parent);
177 }
while (!update_cone());
180 return proceed_with_stability();
198 int Cstable_cones::init_cone(){
206 prepare_cocircular_lists();
214 centre = vicinity[first_cone];
216 centre_idx = first_cone;
221 compute_cone_contents();
233 int Cstable_cones::test_cone(){
234 Creference weighted_cone_ref;
251 cone_candidate = cone;
252 if (cone.ref.not_empty()){
253 hc->insert(&cone_candidate, parent, child,
false,
false);
256 cone_candidate = cone;
257 cone_candidate+= *parent + *child;
258 hc->insert(&cone_candidate, parent, child,
true,
true);
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(){
282 if (centre_idx==vicinity_size)
284 if (centre_idx==first_cone)
296 centre->is_inside->cone =
true;
299 dpt += fabs(child->px)+fabs(child->py);
303 centre = vicinity[centre_idx];
311 if (cocircular_check())
312 return update_cone();
319 if ((centre->side) && (cone.ref.not_empty())){
324 centre->is_inside->cone =
false;
327 dpt += fabs(child->px)+fabs(child->py);
335 if ((dpt>
PT_TSHOLD*(fabs(cone.px)+fabs(cone.py))) && (cone.ref.not_empty())){
336 recompute_cone_contents();
338 if (cone.ref.is_empty()){
352 int Cstable_cones::proceed_with_stability(){
357 for (i=0;i<=hc->mask;i++){
359 elm = hc->hash_array[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
391 nb_hash_cones = hc->n_cones;
392 nb_hash_occupied = hc->n_occupied_cells;
398 return protocones.size();
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)) {
543 if (stable) hc->insert(&candidate);
553 bool Cstable_cones::cocircular_check(){
561 if (centre->cocircular.empty())
return false;
564 if ((centre->side) && (cone.ref.not_empty())){
569 centre->is_inside->cone =
false;
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);
588 centre->cocircular.push_back(centre);
592 for(list<Cvicinity_elm *>::iterator it = centre->cocircular.begin();
593 it != centre->cocircular.end(); it++) {
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 >
681 start(vicinity.begin()+first_cone, vicinity.begin(), vicinity.end());
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(){
729 for (i=0;i<vicinity_size;i++){
731 if ((vicinity[i]->side) && (vicinity[i]->is_inside->cone))
732 cone += *vicinity[i]->v;
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))) {
750 if (cone.ref.is_empty()) {
751 this_cone = Cmomentum();
754 this_cone = Cmomentum();
761 for (
unsigned int i=0;i<vicinity_size;i++){
763 if ((vicinity[i]->side) && (vicinity[i]->is_inside->cone))
764 this_cone += *vicinity[i]->v;
791 Creference Cstable_cones::circle_intersect(
double cx,
double cy){
792 Creference intersection;
796 for (i=0;i<n_part;i++){
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);