SISCone  3.0.5
vicinity.cpp
1 // File: vicinity.cpp //
3 // Description: source file for particle vicinity (Cvicinity class) //
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:: 388 $//
24 // $Date:: 2016-03-03 10:42:25 +0100 (Thu, 03 Mar 2016) $//
26 
27 #include "vicinity.h"
28 #include <math.h>
29 #include <algorithm>
30 #include <iostream>
31 
32 namespace siscone{
33 
34 using namespace std;
35 
36 /*************************************************************
37  * Cvicinity_elm implementation *
38  * element in the vicinity of a parent. *
39  * class used to manage one points in the vicinity *
40  * of a parent point. *
41  *************************************************************/
42 
43 // ordering pointers to Cvicinity_elm
44 //------------------------------------
45 bool ve_less(Cvicinity_elm *ve1, Cvicinity_elm *ve2){
46  return ve1->angle < ve2->angle;
47 }
48 
49 
50 /*************************************************************
51  * Cvicinity implementation *
52  * list of element in the vicinity of a parent. *
53  * class used to manage the points which are in the vicinity *
54  * of a parent point. The construction of the list can be *
55  * made from a list of points or from a quadtree. *
56  *************************************************************/
57 
58 // default constructor
59 //---------------------
61  n_part = 0;
62 
63  ve_list = NULL;
64 #ifdef USE_QUADTREE_FOR_STABILITY_TEST
65  quadtree = NULL;
66 #endif
67 
68  parent = NULL;
69  VR2 = VR = 0.0;
70 
71 }
72 
73 // constructor with initialisation
74 //---------------------------------
75 Cvicinity::Cvicinity(vector<Cmomentum> &_particle_list){
76  parent = NULL;
77 #ifdef USE_QUADTREE_FOR_STABILITY_TEST
78  quadtree = NULL;
79 #endif
80  VR2 = VR = 0.0;
81 
82  ve_list = NULL;
83  set_particle_list(_particle_list);
84 }
85 
86 // default destructor
87 //--------------------
89  if (ve_list!=NULL)
90  delete[] ve_list;
91 
92 #ifdef USE_QUADTREE_FOR_STABILITY_TEST
93  if (quadtree!=NULL)
94  delete quadtree;
95 #endif
96 }
97 
98 /*
99  * set the particle_list
100  * - particle_list list of particles (type Cmomentum)
101  * - n number of particles in the list
102  ************************************************************/
103 void Cvicinity::set_particle_list(vector<Cmomentum> &_particle_list){
104  int i,j;
105 #ifdef USE_QUADTREE_FOR_STABILITY_TEST
106  double eta_max=0.0;
107 #endif
108 
109  // if the particle list is not empty, destroy it !
110  if (ve_list!=NULL){
111  delete[] ve_list;
112  }
113  vicinity.clear();
114 #ifdef USE_QUADTREE_FOR_STABILITY_TEST
115  if (quadtree!=NULL)
116  delete quadtree;
117 #endif
118 
119  // allocate memory array for particles
120  // Note: - we compute max for |eta|
121  // - we allocate indices to particles
122  n_part = 0;
123  plist.clear();
124  pincluded.clear();
125  for (i=0;i<(int) _particle_list.size();i++){
126  // if a particle is colinear with the beam (infinite rapidity)
127  // we do not take it into account
128  if (fabs(_particle_list[i].pz)!=_particle_list[i].E){
129  plist.push_back(_particle_list[i]);
130  pincluded.push_back(Cvicinity_inclusion()); // zero inclusion status
131 
132  // the parent_index is handled in the split_merge because
133  // of our multiple-pass procedure.
134  // Hence, it is not required here any longer.
135  // plist[n_part].parent_index = i;
136  plist[n_part].index = n_part;
137 
138  // make sure the reference is randomly created
139  plist[n_part].ref.randomize();
140 
141 #ifdef USE_QUADTREE_FOR_STABILITY_TEST
142  if (fabs(plist[n_part].eta)>eta_max) eta_max=fabs(plist[n_part].eta);
143 #endif
144 
145  n_part++;
146  }
147  }
148 
149  // allocate quadtree and vicinity_elm list
150  // note: we set phi in [-pi:pi] as it is the natural range for atan2!
151  ve_list = new Cvicinity_elm[2*n_part];
152 #ifdef USE_QUADTREE_FOR_STABILITY_TEST
153  eta_max+=0.1;
154  quadtree = new Cquadtree(0.0, 0.0, eta_max, M_PI);
155 #endif
156 
157  // append particle to the vicinity_elm list
158  j = 0;
159  for (i=0;i<n_part;i++){
160 #ifdef USE_QUADTREE_FOR_STABILITY_TEST
161  quadtree->add(&plist[i]);
162 #endif
163  ve_list[j].v = ve_list[j+1].v = &plist[i];
164  ve_list[j].is_inside = ve_list[j+1].is_inside = &(pincluded[i]);
165  j+=2;
166  }
167 
168 }
169 
170 
171 /*
172  * build the vicinity list from a list of points.
173  * - _parent reference particle
174  * - _VR vicinity radius
175  ************************************************************/
176 void Cvicinity::build(Cmomentum *_parent, double _VR){
177  int i;
178 
179  // set parent and radius
180  parent = _parent;
181  VR = _VR;
182  VR2 = VR*VR;
183  R2 = 0.25*VR2;
184  R = 0.5*VR;
185  inv_R_EPS_COCIRC = 1.0 / R / EPSILON_COCIRCULAR;
186  inv_R_2EPS_COCIRC = 0.5 / R / EPSILON_COCIRCULAR;
187 
188  // clear vicinity
189  vicinity.clear();
190 
191  // init parent variables
192  pcx = parent->eta;
193  pcy = parent->phi;
194 
195  // really browse the particle list
196  for (i=0;i<n_part;i++){
197  append_to_vicinity(&plist[i]);
198  }
199 
200  // sort the vicinity
201  sort(vicinity.begin(), vicinity.end(), ve_less);
202 
203  vicinity_size = vicinity.size();
204 }
205 
206 
208 inline double sort_angle(double s, double c){
209  if (s==0) return (c>0) ? 0.0 : 2.0;
210  double t=c/s;
211  return (s>0) ? 1-t/(1+fabs(t)) : 3-t/(1+fabs(t));
212 }
213 
214 
215 /*
216  * append a particle to the 'vicinity' list after
217  * having computed the angular-ordering quantities
218  * - v vector to test
219  **********************************************************/
221  double dx, dy, d2;
222 
223  // skip the particle itself)
224  if (v==parent)
225  return;
226 
227  int i=2*(v->index);
228 
229  // compute the distance of the i-th particle with the parent
230  dx = v->eta - pcx;
231  dy = v->phi - pcy;
232 
233  // pay attention to the periodicity in phi !
234  if (dy>M_PI)
235  dy -= twopi;
236  else if (dy<-M_PI)
237  dy += twopi;
238 
239  d2 = dx*dx+dy*dy;
240 
241  // really check if the distance is less than VR
242  if (d2<VR2){
243  double s,c,tmp;
244 
245  // compute the angles used for future ordering ...
246  // - build temporary variables used for the computation
247  //d = sqrt(d2);
248  tmp = sqrt(VR2/d2-1);
249 
250  // first angle (+)
251  c = 0.5*(dx-dy*tmp); // cosine of (parent,child) pair w.r.t. horizontal
252  s = 0.5*(dy+dx*tmp); // sine of (parent,child) pair w.r.t. horizontal
253  ve_list[i].angle = sort_angle(s,c);
254  ve_list[i].eta = pcx+c;
255  ve_list[i].phi = phi_in_range(pcy+s);
256  ve_list[i].side = true;
257  ve_list[i].cocircular.clear();
258  vicinity.push_back(&(ve_list[i]));
259 
260  // second angle (-)
261  c = 0.5*(dx+dy*tmp); // cosine of (parent,child) pair w.r.t. horizontal
262  s = 0.5*(dy-dx*tmp); // sine of (parent,child) pair w.r.t. horizontal
263  ve_list[i+1].angle = sort_angle(s,c);
264  ve_list[i+1].eta = pcx+c;
265  ve_list[i+1].phi = phi_in_range(pcy+s);
266  ve_list[i+1].side = false;
267  ve_list[i+1].cocircular.clear();
268  vicinity.push_back(&(ve_list[i+1]));
269 
270  // now work out the cocircularity range for the two points (range
271  // of angle within which the points stay within a distance
272  // EPSILON_COCIRCULAR of circule
273  // P = parent; C = child; O = Origin (center of circle)
274  Ctwovect OP(pcx - ve_list[i+1].eta, phi_in_range(pcy-ve_list[i+1].phi));
275  Ctwovect OC(v->eta - ve_list[i+1].eta,
276  phi_in_range(v->phi-ve_list[i+1].phi));
277 
278  // two sources of error are (GPS CCN29-19) epsilon/(R sin theta)
279  // and sqrt(2*epsilon/(R (1-cos theta))) and the way things work
280  // out, it is the _smaller_ of the two that is relevant [NB have
281  // changed definition of theta here relative to that used in
282  // CCN29] [NB2: write things so as to avoid zero denominators and
283  // to minimize the multiplications, divisions and above all sqrts
284  // -- that means that c & s are defined including a factor of VR2]
285  c = dot_product(OP,OC);
286  s = fabs(cross_product(OP,OC));
287  double inv_err1 = s * inv_R_EPS_COCIRC;
288  double inv_err2_sq = (R2-c) * inv_R_2EPS_COCIRC;
289  ve_list[i].cocircular_range = pow2(inv_err1) > inv_err2_sq ?
290  1.0/inv_err1 :
291  sqrt(1.0/inv_err2_sq);
292  ve_list[i+1].cocircular_range = ve_list[i].cocircular_range;
293  }
294 }
295 
296 }
base class for dynamic coordinates management
Definition: momentum.h:49
int index
internal particle number
Definition: momentum.h:117
double eta
particle pseudo-rapidity
Definition: momentum.h:114
double phi
particle azimuthal angle
Definition: momentum.h:115
Implementation of a 2D quadtree.
Definition: quadtree.h:43
class for holding a two-vector
Definition: geom_2d.h:73
element in the vicinity of a parent.
Definition: vicinity.h:63
a class to keep track of inclusion status in cone and in cocircular region while using minimal resour...
Definition: vicinity.h:46
void append_to_vicinity(Cmomentum *v)
append a particle to the 'vicinity' list after having tested it and computed the angular-ordering qua...
Definition: vicinity.cpp:220
void build(Cmomentum *_parent, double _VR)
build the vicinity list from the list of points.
Definition: vicinity.cpp:176
Cvicinity()
default constructor
Definition: vicinity.cpp:60
~Cvicinity()
default destructor
Definition: vicinity.cpp:88
void set_particle_list(std::vector< Cmomentum > &_particle_list)
set the particle_list
Definition: vicinity.cpp:103
#define EPSILON_COCIRCULAR
The following parameter controls cocircular situations.
Definition: defines.h:82
const double twopi
definition of 2*M_PI which is useful a bit everyhere!
Definition: defines.h:114