SISCone  2.0.5
protocones.cpp
1 
2 // File: protocones.cpp //
3 // Description: source file for stable cones determination (Cstable_cones) //
4 // This file is part of the SISCone project. //
5 // For more details, see http://projects.hepforge.org/siscone //
6 // //
7 // Copyright (c) 2006 Gavin Salam and Gregory Soyez //
8 // //
9 // This program is free software; you can redistribute it and/or modify //
10 // it under the terms of the GNU General Public License as published by //
11 // the Free Software Foundation; either version 2 of the License, or //
12 // (at your option) any later version. //
13 // //
14 // This program is distributed in the hope that it will be useful, //
15 // but WITHOUT ANY WARRANTY; without even the implied warranty of //
16 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the //
17 // GNU General Public License for more details. //
18 // //
19 // You should have received a copy of the GNU General Public License //
20 // along with this program; if not, write to the Free Software //
21 // Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301, USA //
22 // //
23 // $Revision:: 322 $//
24 // $Date:: 2011-11-15 10:12:36 +0100 (Tue, 15 Nov 2011) $//
26 
27 /*******************************************************
28  * Introductory note: *
29  * Since this file has many member functions, we have *
30  * structured them in categories: *
31  * INITIALISATION FUNCTIONS *
32  * - ctor() *
33  * - ctor(particle_list) *
34  * - dtor() *
35  * - init(particle_list) *
36  * ALGORITHM MAIN ENTRY *
37  * - get_stable_cone(radius) *
38  * ALGORITHM MAIN STEPS *
39  * - init_cone() *
40  * - test_cone() *
41  * - update_cone() *
42  * - proceed_with_stability() *
43  * ALGORITHM MAIN STEPS FOR COCIRCULAR SITUATIONS *
44  * - cocircular_pt_less(v1, v2) *
45  * - prepare_cocircular_list() *
46  * - test_cone_cocircular() *
47  * - test_stability(candidate, border_list) *
48  * - updat_cone_cocircular() *
49  * RECOMPUTATION OF CONE CONTENTS *
50  * - compute_cone_contents() *
51  * - recompute_cone_contents() *
52  * - recompute_cone_contents_if_needed() *
53  * VARIOUS TOOLS *
54  * - circle_intersect() *
55  * - is_inside() *
56  * - abs_dangle() *
57  *******************************************************/
58 
59 #include "protocones.h"
60 #include "siscone_error.h"
61 #include "defines.h"
62 #include <math.h>
63 #include <iostream>
64 #include "circulator.h"
65 #include <algorithm>
66 
67 namespace siscone{
68 
69 using namespace std;
70 
71 /**********************************************************************
72  * Cstable_cones implementation *
73  * Computes the list of stable comes from a particle list. *
74  * This class does the first fundamental task of te cone algorithm: *
75  * it is used to compute the list of stable cones given a list *
76  * of particles. *
77  **********************************************************************/
78 
80 // INITIALISATION FUNCTIONS //
81 // - ctor() //
82 // - ctor(particle_list) //
83 // - dtor() //
84 // - init(particle_list) //
86 
87 // default ctor
88 //--------------
90  nb_tot = 0;
91  hc = NULL;
92 }
93 
94 // ctor with initialisation
95 //--------------------------
96 Cstable_cones::Cstable_cones(vector<Cmomentum> &_particle_list)
97  : Cvicinity(_particle_list){
98 
99  nb_tot = 0;
100  hc = NULL;
101 }
102 
103 // default dtor
104 //--------------
106  if (hc!=NULL) delete hc;
107 }
108 
109 /*
110  * initialisation
111  * - _particle_list list of particles
112  * - _n number of particles
113  *********************************************************************/
114 void Cstable_cones::init(vector<Cmomentum> &_particle_list){
115  // check already allocated mem
116  if (hc!=NULL){
117  delete hc;
118  }
119  if (protocones.size()!=0)
120  protocones.clear();
121 
122  multiple_centre_done.clear();
123 
124  // initialisation
125  set_particle_list(_particle_list);
126 }
127 
128 
130 // ALGORITHM MAIN ENTRY //
131 // - get_stable_cone(radius) //
133 
134 /*
135  * compute stable cones.
136  * This function really does the job i.e. computes
137  * the list of stable cones (in a seedless way)
138  * - _radius: radius of the cones
139  * The number of stable cones found is returned
140  *********************************************************************/
142  int p_idx;
143 
144  // check if everything is correctly initialised
145  if (n_part==0){
146  return 0;
147  }
148 
149  R = _radius;
150  R2 = R*R;
151 
152  // allow hash for cones candidates
153  hc = new hash_cones(n_part, R2);
154 
155  // browse all particles
156  for (p_idx=0;p_idx<n_part;p_idx++){
157  // step 0: compute the child list CL.
158  // Note that this automatically sets the parent P
159  build(&plist[p_idx], 2.0*R);
160 
161  // special case:
162  // if the vicinity is empty, the parent particle is a
163  // stable cone by itself. Add it to protocones list.
164  if (vicinity_size==0){
165  protocones.push_back(*parent);
166  continue;
167  }
168 
169  // step 1: initialise with the first cone candidate
170  init_cone();
171 
172  do{
173  // step 2: test cone stability for that pair (P,C)
174  test_cone();
175 
176  // step 3: go to the next cone child candidate C
177  } while (!update_cone());
178  }
179 
180  return proceed_with_stability();
181 }
182 
183 
185 // ALGORITHM MAIN STEPS //
186 // - init_cone() //
187 // - test_cone() //
188 // - update_cone() //
189 // - proceed_with_stability() //
191 
192 /*
193  * initialise the cone.
194  * We take the first particle in the angular ordering to compute
195  * this one
196  * return 0 on success, 1 on error
197  *********************************************************************/
198 int Cstable_cones::init_cone(){
199  // The previous version of the algorithm was starting the
200  // loop around vicinity elements with the "most isolated" child.
201  // given the nodist method to calculate the cone contents, we no
202  // longer need to worry about which cone comes first...
203  first_cone=0;
204 
205  // now make sure we have lists of the cocircular particles
206  prepare_cocircular_lists();
207 
208  //TODO? deal with a configuration with only degeneracies ?
209  // The only possibility seems a regular hexagon with a parent point
210  // in the centre. And this situation is by itself unclear.
211  // Hence, we do nothing here !
212 
213  // init set child C
214  centre = vicinity[first_cone];
215  child = centre->v;
216  centre_idx = first_cone;
217 
218  // build the initial cone (nodist: avoids calculating distances --
219  // just deduces contents by circulating around all in/out operations)
220  // this function also sets the list of included particles
221  compute_cone_contents();
222 
223  return 0;
224 }
225 
226 
227 /*
228  * test cones.
229  * We check if the cone(s) built with the present parent and child
230  * are stable
231  * return 0 on success 1 on error
232  *********************************************************************/
233 int Cstable_cones::test_cone(){
234  Creference weighted_cone_ref;
235 
236  // depending on the side we are taking the child particle,
237  // we test different configuration.
238  // Each time, two configurations are tested in such a way that
239  // all 4 possible cases (parent or child in or out the cone)
240  // are tested when taking the pair of particle parent+child
241  // and child+parent.
242 
243  // here are the tests entering the first series:
244  // 1. check if the cone is already inserted
245  // 2. check cone stability for the parent and child particles
246 
247  if (centre->side){
248  // test when both particles are not in the cone
249  // or when both are in.
250  // Note: for the totally exclusive case, test emptyness before
251  cone_candidate = cone;
252  if (cone.ref.not_empty()){
253  hc->insert(&cone_candidate, parent, child, false, false);
254  }
255 
256  cone_candidate = cone;
257  cone_candidate+= *parent + *child;
258  hc->insert(&cone_candidate, parent, child, true, true);
259  } else {
260  // test when 1! of the particles is in the cone
261  cone_candidate = cone + *parent;
262  hc->insert(&cone_candidate, parent, child, true, false);
263 
264  cone_candidate = cone + *child;
265  hc->insert(&cone_candidate, parent, child, false, true);
266  }
267 
268  nb_tot+=2;
269 
270  return 0;
271 }
272 
273 
274 /*
275  * update the cone
276  * go to the next child for that parent and update 'cone' appropriately
277  * return 0 if update candidate found, 1 otherwise
278  ***********************************************************************/
279 int Cstable_cones::update_cone(){
280  // get the next child and centre
281  centre_idx++;
282  if (centre_idx==vicinity_size)
283  centre_idx=0;
284  if (centre_idx==first_cone)
285  return 1;
286 
287  // update the cone w.r.t. the old child
288  // only required if the old child is entering inside in which
289  // case we need to add it. We also know that the child is
290  // inside iff its side is -.
291  if (!centre->side){
292  // update cone
293  cone += (*child);
294 
295  // update info on particles inside
296  centre->is_inside->cone = true;
297 
298  // update stability check quantities
299  dpt += fabs(child->px)+fabs(child->py);
300  }
301 
302  // update centre and child to correspond to the new position
303  centre = vicinity[centre_idx];
304  child = centre->v;
305 
306  // check cocircularity
307  // note that if cocirculaity is detected (i.e. if we receive 1
308  // in the next test), we need to recall 'update_cone' directly
309  // since tests and remaining part of te update has been performed
310  //if (cocircular_check())
311  if (cocircular_check())
312  return update_cone();
313 
314 
315  // update the cone w.r.t. the new child
316  // only required if the new child was already inside in which
317  // case we need to remove it. We also know that the child is
318  // inside iff its side is +.
319  if ((centre->side) && (cone.ref.not_empty())){
320  // update cone
321  cone -= (*child);
322 
323  // update info on particles inside
324  centre->is_inside->cone = false;
325 
326  // update stability check quantities
327  dpt += fabs(child->px)+fabs(child->py); //child->perp2();
328  }
329 
330  // check that the addition and subtraction of vectors does
331  // not lead to too much rounding error
332  // for that, we compute the sum of pt modifications and of |pt|
333  // since last recomputation and once the ratio overpasses a threshold
334  // we recompute vicinity.
335  if ((dpt>PT_TSHOLD*(fabs(cone.px)+fabs(cone.py))) && (cone.ref.not_empty())){
336  recompute_cone_contents();
337  }
338  if (cone.ref.is_empty()){
339  cone = Cmomentum();
340  dpt=0.0;
341  }
342 
343  return 0;
344 }
345 
346 
347 /*
348  * compute stability of all enumerated candidates.
349  * For all candidate cones which are stable w.r.t. their border particles,
350  * pass the last test: stability with quadtree intersection
351  ************************************************************************/
352 int Cstable_cones::proceed_with_stability(){
353  int i; // ,n;
354  hash_element *elm;
355 
356  //n=0;
357  for (i=0;i<=hc->mask;i++){
358  // test ith cell of the hash array
359  elm = hc->hash_array[i];
360 
361  // browse elements therein
362  while (elm!=NULL){
363  // test stability
364  if (elm->is_stable){
365  // stability is not ensured by all pairs of "edges" already browsed
366 #ifdef USE_QUADTREE_FOR_STABILITY_TEST
367  // => testing stability with quadtree intersection
368  if (quadtree->circle_intersect(elm->eta, elm->phi, R2)==elm->ref){
369 #else
370  // => testing stability with the particle-list intersection
371  if (circle_intersect(elm->eta, elm->phi)==elm->ref){
372 #endif
373  // add it to the list of protocones
374  // note that in its present form, we do not allocate the
375  // 4-vector components of the momentum. There's no need to
376  // do it here as it will be recomputed in
377  // Csplit_merge::add_protocones
378  protocones.push_back(Cmomentum(elm->eta, elm->phi, elm->ref));
379  }
380  }
381 
382  // jump to the next one
383  elm = elm->next;
384  }
385  }
386 
387  // free hash
388  // we do that at this level because hash eats rather a lot of memory
389  // we want to free it before running the split/merge algorithm
390 #ifdef DEBUG_STABLE_CONES
391  nb_hash_cones = hc->n_cones;
392  nb_hash_occupied = hc->n_occupied_cells;
393 #endif
394 
395  delete hc;
396  hc=NULL;
397 
398  return protocones.size();
399 }
400 
401 
403 // ALGORITHM MAIN STEPS FOR COCIRCULAR SITUATIONS //
404 // - cocircular_pt_less(v1, v2) //
405 // - prepare_cocircular_list() //
406 // - test_cone_cocircular() //
407 // - test_stability(candidate, border_vect) //
408 // - updat_cone_cocircular() //
410 
412 bool cocircular_pt_less(Cmomentum *v1, Cmomentum *v2){
413  return v1->perp2() < v2->perp2();
414 }
415 
416 /*
417  * run through the vicinity of the current parent and for each child
418  * establish which other members are cocircular... Note that the list
419  * associated with each child contains references to vicinity
420  * elements: thus two vicinity elements each associated with one given
421  * particle may appear in a list -- this needs to be watched out for
422  * later on...
423  **********************************************************************/
424 void Cstable_cones::prepare_cocircular_lists() {
425  circulator<vector<Cvicinity_elm*>::iterator > here(vicinity.begin(),
426  vicinity.begin(),
427  vicinity.end());
428 
429  circulator<vector<Cvicinity_elm*>::iterator > search(here);
430 
431  do {
432  Cvicinity_elm* here_pntr = *here();
433  search.set_position(here);
434 
435  // search forwards for things that should have "here" included in
436  // their cocircularity list
437  while (true) {
438  ++search;
439  if ( abs_dphi((*search())->angle, here_pntr->angle) <
440  here_pntr->cocircular_range
441  && search() != here()) {
442  (*search())->cocircular.push_back(here_pntr);
443  } else {
444  break;
445  }
446  }
447 
448  // search backwards
449  search.set_position(here);
450  while (true) {
451  --search;
452  if ( abs_dphi((*search())->angle, here_pntr->angle) <
453  here_pntr->cocircular_range
454  && search() != here()) {
455  (*search())->cocircular.push_back(here_pntr);
456  } else {
457  break;
458  }
459  }
460 
461  ++here;
462  } while (here() != vicinity.begin());
463 
464 }
465 
466 /*
467  * Testing cocircular configurations in p^3 time,
468  * rather than 2^p time; we will test all contiguous subsets of points
469  * on the border --- note that this is till probably overkill, since
470  * in principle we only have to test situations where up to a
471  * half-circle is filled (but going to a full circle is simpler)
472  ******************************************************************/
473 void Cstable_cones::test_cone_cocircular(Cmomentum & borderless_cone,
474  list<Cmomentum *> & border_list) {
475  vector<Cborder_store> border_vect;
476 
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));
481  }
482 
483  // get them into order of angle
484  sort(border_vect.begin(), border_vect.end());
485 
486  // set up some circulators, since these will help us go around the
487  // circle easily
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);
491 
492  // test the borderless cone
493  Cmomentum candidate = borderless_cone;
494  candidate.build_etaphi();
495  if (candidate.ref.not_empty())
496  test_stability(candidate, border_vect);
497 
498  do {
499  // reset status wrt inclusion in the cone
500  mid = start;
501  do {
502  mid()->is_in = false;
503  } while (++mid != start);
504 
505  // now run over all inclusion possibilities with this starting point
506  candidate = borderless_cone;
507  while (++mid != start) {
508  // will begin with start+1 and go up to start-1
509  mid()->is_in = true;
510  candidate += *(mid()->mom);
511  test_stability(candidate, border_vect);
512  }
513 
514  } while (++start != end);
515 
516  // mid corresponds to momentum that we need to include to get the
517  // full cone
518  mid()->is_in = true;
519  candidate += *(mid()->mom);
520  test_stability(candidate, border_vect);
521 }
522 
523 
530 void Cstable_cones::test_stability(Cmomentum & candidate, const vector<Cborder_store> & border_vect) {
531 
532  // this almost certainly has not been done...
533  candidate.build_etaphi();
534 
535  bool stable = true;
536  for (unsigned i = 0; i < border_vect.size(); i++) {
537  if (is_inside(&candidate, border_vect[i].mom) ^ (border_vect[i].is_in)) {
538  stable = false;
539  break; // it's unstable so there's no point continuing
540  }
541  }
542 
543  if (stable) hc->insert(&candidate);
544 }
545 
546 /*
547  * check if we are in a situation of cocircularity.
548  * if it is the case, update and test in the corresponding way
549  * return 'false' if no cocircularity detected, 'true' otherwise
550  * Note that if cocircularity is detected, we need to
551  * recall 'update' from 'update' !!!
552  ***************************************************************/
553 bool Cstable_cones::cocircular_check(){
554  // check if many configurations have the same centre.
555  // if this is the case, branch on the algorithm for this
556  // special case.
557  // Note that those situation, being considered separately in
558  // test_cone_multiple, must only be considered here if all
559  // angles are on the same side (this avoid multiple counting)
560 
561  if (centre->cocircular.empty()) return false;
562 
563  // first get cone into status required at end...
564  if ((centre->side) && (cone.ref.not_empty())){
565  // update cone
566  cone -= (*child);
567 
568  // update info on particles inside
569  centre->is_inside->cone = false;
570 
571  // update stability check quantities
572  dpt += fabs(child->px)+fabs(child->py); //child->perp2();
573  }
574 
575 
576  // now establish the list of unique children in the list
577  // first make sure parent and child are in!
578 
579  list<Cvicinity_inclusion *> removed_from_cone;
580  list<Cvicinity_inclusion *> put_in_border;
581  list<Cmomentum *> border_list;
582 
583  Cmomentum cone_removal;
584  Cmomentum border = *parent;
585  border_list.push_back(parent);
586 
587  // make sure child appears in the border region
588  centre->cocircular.push_back(centre);
589 
590  // now establish the full contents of the cone minus the cocircular
591  // region and of the cocircular region itself
592  for(list<Cvicinity_elm *>::iterator it = centre->cocircular.begin();
593  it != centre->cocircular.end(); it++) {
594 
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);
599  }
600 
601  // if a point appears twice (i.e. with + and - sign) in the list of
602  // points on the border, we take care not to include it twice.
603  // Note that this situation may appear when a point is at a distance
604  // close to 2R from the parent
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);
610  }
611  }
612 
613 
614  // figure out whether this pairing has been observed before
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))
621  consider = false;
622  }
623 
624  // now prepare the hard work
625  if (consider) {
626  // record the fact that we've now seen this combination
627  multiple_centre_done.push_back(pair<Creference,Creference>(borderless_cone.ref,
628  border.ref));
629 
630  // first figure out whether our cone momentum is good
631  double local_dpt = fabs(cone_removal.px) + fabs(cone_removal.py);
632  double total_dpt = dpt + local_dpt;
633 
634  recompute_cone_contents_if_needed(borderless_cone, total_dpt);
635  if (total_dpt == 0) {
636  // a recomputation has taken place -- so take advantage of this
637  // and update the member cone momentum
638  cone = borderless_cone + cone_removal;
639  dpt = local_dpt;
640  }
641 
642  test_cone_cocircular(borderless_cone, border_list);
643  }
644 
645 
646  // relabel things that were in the cone but got removed
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;
650  }
651 
652  // relabel things that got put into the border
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;
656  }
657 
658  // we're done with everything -- return true to signal to user that we've
659  // been through the co-circularity rigmarole
660  return true;
661 }
662 
663 
665 // RECOMPUTATION OF CONE CONTENTS //
666 // - compute_cone_contents() //
667 // - recompute_cone_contents() //
668 // - recompute_cone_contents_if_needed() //
670 
679 void Cstable_cones::compute_cone_contents() {
680  circulator<vector<Cvicinity_elm*>::iterator >
681  start(vicinity.begin()+first_cone, vicinity.begin(), vicinity.end());
682 
683  circulator<vector<Cvicinity_elm*>::iterator > here(start);
684 
685  // note that in the following algorithm, the cone contents never includes
686  // the child. Indeed, if it has positive sign, then it will be set as
687  // outside at the last step in the loop. If it has negative sign, then the
688  // loop will at some point go to the corresponding situation with positive
689  // sign and set the inclusion status to 0.
690 
691  do {
692  // as we leave this position a particle enters if its side is
693  // negative (i.e. the centre is the one at -ve angle wrt to the
694  // parent-child line
695  if (!(*here())->side) ((*here())->is_inside->cone) = 1;
696 
697  // move on to the next position
698  ++here;
699 
700  // as we arrive at this position a particle leaves if its side is positive
701  if ((*here())->side) ((*here())->is_inside->cone) = 0;
702  } while (here != start);
703 
704  // once we've reached the start the 'is_inside' information should be
705  // 100% complete, so we can use it to calculate the cone contents
706  // and then exit
707  recompute_cone_contents();
708  return;
709 
710 }
711 
712 
713 /*
714  * compute the cone momentum from particle list.
715  * in this version, we use the 'pincluded' information
716  * from the Cvicinity class
717  */
718 void Cstable_cones::recompute_cone_contents(){
719  unsigned int i;
720 
721  // set momentum to 0
722  cone = Cmomentum();
723 
724  // Important note: we can browse only the particles
725  // in vicinity since all particles in the cone are
726  // withing a distance 2R w.r.t. parent hence in vicinity.
727  // Among those, we only add the particles for which 'is_inside' is true !
728  // This methos rather than a direct comparison avoids rounding errors
729  for (i=0;i<vicinity_size;i++){
730  // to avoid double-counting, only use particles with + angle
731  if ((vicinity[i]->side) && (vicinity[i]->is_inside->cone))
732  cone += *vicinity[i]->v;
733  }
734 
735  // set check variables back to 0
736  dpt = 0.0;
737 }
738 
739 
740 /*
741  * if we have gone beyond the acceptable threshold of change, compute
742  * the cone momentum from particle list. in this version, we use the
743  * 'pincluded' information from the Cvicinity class, but we don't
744  * change the member cone, only the locally supplied one
745  */
746 void Cstable_cones::recompute_cone_contents_if_needed(Cmomentum & this_cone,
747  double & this_dpt){
748 
749  if (this_dpt > PT_TSHOLD*(fabs(this_cone.px)+fabs(this_cone.py))) {
750  if (cone.ref.is_empty()) {
751  this_cone = Cmomentum();
752  } else {
753  // set momentum to 0
754  this_cone = Cmomentum();
755 
756  // Important note: we can browse only the particles
757  // in vicinity since all particles in the this_cone are
758  // withing a distance 2R w.r.t. parent hence in vicinity.
759  // Among those, we only add the particles for which 'is_inside' is true !
760  // This methos rather than a direct comparison avoids rounding errors
761  for (unsigned int i=0;i<vicinity_size;i++){
762  // to avoid double-counting, only use particles with + angle
763  if ((vicinity[i]->side) && (vicinity[i]->is_inside->cone))
764  this_cone += *vicinity[i]->v;
765  }
766 
767  }
768  // set check variables back to 0
769  this_dpt = 0.0;
770  }
771 
772 }
773 
774 
776 // VARIOUS TOOLS //
777 // - circle_intersect() //
778 // - is_inside() //
779 // - abs_dangle() //
781 
782 
783 /*
784  * circle intersection.
785  * computes the intersection with a circle of given centre and radius.
786  * The output takes the form of a checkxor of the intersection's particles
787  * - cx circle centre x coordinate
788  * - cy circle centre y coordinate
789  * return the checkxor for the intersection
790  ******************************************************************/
791 Creference Cstable_cones::circle_intersect(double cx, double cy){
792  Creference intersection;
793  int i;
794  double dx, dy;
795 
796  for (i=0;i<n_part;i++){
797  // compute the distance of the i-th particle with the parent
798  dx = plist[i].eta - cx;
799  dy = fabs(plist[i].phi - cy);
800 
801  // pay attention to the periodicity in phi !
802  if (dy>M_PI)
803  dy -= twopi;
804 
805  // really check if the distance is less than VR
806  if (dx*dx+dy*dy<R2)
807  intersection+=plist[i].ref;
808  }
809 
810  return intersection;
811 }
812 
813 /*
814  * test if a particle is inside a cone of given centre.
815  * check if the particle of coordinates 'v' is inside the circle of radius R
816  * centered at 'centre'.
817  * - centre centre of the circle
818  * - v particle to test
819  * return true if inside, false if outside
820  *****************************************************************************/
821 inline bool Cstable_cones::is_inside(Cmomentum *centre_in, Cmomentum *v){
822  double dx, dy;
823 
824  dx = centre_in->eta - v->eta;
825  dy = fabs(centre_in->phi - v->phi);
826  if (dy>M_PI)
827  dy -= twopi;
828 
829  return dx*dx+dy*dy<R2;
830 }
831 
832 /*
833  * compute the absolute value of the difference between 2 angles.
834  * We take care of the 2pi periodicity
835  * - angle1 first angle
836  * - angle2 second angle
837  * return the absolute value of the difference between the angles
838  *****************************************************************/
839 inline double abs_dangle(double &angle1, double &angle2){
840  double dphi;
841 
842  dphi = fabs(angle1-angle2);
843  if (dphi>M_PI)
844  dphi = dphi-twopi;
845 
846  return dphi;
847 }
848 
849 }
The SISCone project has been developed by Gavin Salam and Gregory Soyez
Documentation generated on Wed Mar 19 2014 21:30:25 for SISCone by  Doxygen 1.8.1.2