SISCone  2.0.6
vicinity.cpp
1 
2 // 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:: 123 $//
24 // $Date:: 2007-03-01 02:52:16 +0100 (Thu, 01 Mar 2007) $//
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  set_particle_list(_particle_list);
83 }
84 
85 // default destructor
86 //--------------------
88  if (ve_list!=NULL)
89  delete[] ve_list;
90 
91 #ifdef USE_QUADTREE_FOR_STABILITY_TEST
92  if (quadtree!=NULL)
93  delete quadtree;
94 #endif
95 }
96 
97 /*
98  * set the particle_list
99  * - particle_list list of particles (type Cmomentum)
100  * - n number of particles in the list
101  ************************************************************/
102 void Cvicinity::set_particle_list(vector<Cmomentum> &_particle_list){
103  int i,j;
104 #ifdef USE_QUADTREE_FOR_STABILITY_TEST
105  double eta_max=0.0;
106 #endif
107 
108  // if the particle list is not empty, destroy it !
109  if (ve_list!=NULL){
110  delete[] ve_list;
111  }
112  vicinity.clear();
113 #ifdef USE_QUADTREE_FOR_STABILITY_TEST
114  if (quadtree!=NULL)
115  delete quadtree;
116 #endif
117 
118  // allocate memory array for particles
119  // Note: - we compute max for |eta|
120  // - we allocate indices to particles
121  n_part = 0;
122  plist.clear();
123  pincluded.clear();
124  for (i=0;i<(int) _particle_list.size();i++){
125  // if a particle is colinear with the beam (infinite rapidity)
126  // we do not take it into account
127  if (fabs(_particle_list[i].pz)!=_particle_list[i].E){
128  plist.push_back(_particle_list[i]);
129  pincluded.push_back(Cvicinity_inclusion()); // zero inclusion status
130 
131  // the parent_index is handled in the split_merge because
132  // of our multiple-pass procedure.
133  // Hence, it is not required here any longer.
134  // plist[n_part].parent_index = i;
135  plist[n_part].index = n_part;
136 
137  // make sure the reference is randomly created
138  plist[n_part].ref.randomize();
139 
140 #ifdef USE_QUADTREE_FOR_STABILITY_TEST
141  if (fabs(plist[n_part].eta)>eta_max) eta_max=fabs(plist[n_part].eta);
142 #endif
143 
144  n_part++;
145  }
146  }
147 
148  // allocate quadtree and vicinity_elm list
149  // note: we set phi in [-pi:pi] as it is the natural range for atan2!
150  ve_list = new Cvicinity_elm[2*n_part];
151 #ifdef USE_QUADTREE_FOR_STABILITY_TEST
152  eta_max+=0.1;
153  quadtree = new Cquadtree(0.0, 0.0, eta_max, M_PI);
154 #endif
155 
156  // append particle to the vicinity_elm list
157  j = 0;
158  for (i=0;i<n_part;i++){
159 #ifdef USE_QUADTREE_FOR_STABILITY_TEST
160  quadtree->add(&plist[i]);
161 #endif
162  ve_list[j].v = ve_list[j+1].v = &plist[i];
163  ve_list[j].is_inside = ve_list[j+1].is_inside = &(pincluded[i]);
164  j+=2;
165  }
166 
167 }
168 
169 
170 /*
171  * build the vicinity list from a list of points.
172  * - _parent reference particle
173  * - _VR vicinity radius
174  ************************************************************/
175 void Cvicinity::build(Cmomentum *_parent, double _VR){
176  int i;
177 
178  // set parent and radius
179  parent = _parent;
180  VR = _VR;
181  VR2 = VR*VR;
182  R2 = 0.25*VR2;
183  R = 0.5*VR;
184  inv_R_EPS_COCIRC = 1.0 / R / EPSILON_COCIRCULAR;
185  inv_R_2EPS_COCIRC = 0.5 / R / EPSILON_COCIRCULAR;
186 
187  // clear vicinity
188  vicinity.clear();
189 
190  // init parent variables
191  pcx = parent->eta;
192  pcy = parent->phi;
193 
194  // really browse the particle list
195  for (i=0;i<n_part;i++){
196  append_to_vicinity(&plist[i]);
197  }
198 
199  // sort the vicinity
200  sort(vicinity.begin(), vicinity.end(), ve_less);
201 
202  vicinity_size = vicinity.size();
203 }
204 
205 
207 inline double sort_angle(double s, double c){
208  if (s==0) return (c>0) ? 0.0 : 2.0;
209  double t=c/s;
210  return (s>0) ? 1-t/(1+fabs(t)) : 3-t/(1+fabs(t));
211 }
212 
213 
214 /*
215  * append a particle to the 'vicinity' list after
216  * having computed the angular-ordering quantities
217  * - v vector to test
218  **********************************************************/
220  double dx, dy, d2;
221 
222  // skip the particle itself)
223  if (v==parent)
224  return;
225 
226  int i=2*(v->index);
227 
228  // compute the distance of the i-th particle with the parent
229  dx = v->eta - pcx;
230  dy = v->phi - pcy;
231 
232  // pay attention to the periodicity in phi !
233  if (dy>M_PI)
234  dy -= twopi;
235  else if (dy<-M_PI)
236  dy += twopi;
237 
238  d2 = dx*dx+dy*dy;
239 
240  // really check if the distance is less than VR
241  if (d2<VR2){
242  double s,c,tmp;
243 
244  // compute the angles used for future ordering ...
245  // - build temporary variables used for the computation
246  //d = sqrt(d2);
247  tmp = sqrt(VR2/d2-1);
248 
249  // first angle (+)
250  c = 0.5*(dx-dy*tmp); // cosine of (parent,child) pair w.r.t. horizontal
251  s = 0.5*(dy+dx*tmp); // sine of (parent,child) pair w.r.t. horizontal
252  ve_list[i].angle = sort_angle(s,c);
253  ve_list[i].eta = pcx+c;
254  ve_list[i].phi = phi_in_range(pcy+s);
255  ve_list[i].side = true;
256  ve_list[i].cocircular.clear();
257  vicinity.push_back(&(ve_list[i]));
258 
259  // second angle (-)
260  c = 0.5*(dx+dy*tmp); // cosine of (parent,child) pair w.r.t. horizontal
261  s = 0.5*(dy-dx*tmp); // sine of (parent,child) pair w.r.t. horizontal
262  ve_list[i+1].angle = sort_angle(s,c);
263  ve_list[i+1].eta = pcx+c;
264  ve_list[i+1].phi = phi_in_range(pcy+s);
265  ve_list[i+1].side = false;
266  ve_list[i+1].cocircular.clear();
267  vicinity.push_back(&(ve_list[i+1]));
268 
269  // now work out the cocircularity range for the two points (range
270  // of angle within which the points stay within a distance
271  // EPSILON_COCIRCULAR of circule
272  // P = parent; C = child; O = Origin (center of circle)
273  Ctwovect OP(pcx - ve_list[i+1].eta, phi_in_range(pcy-ve_list[i+1].phi));
274  Ctwovect OC(v->eta - ve_list[i+1].eta,
275  phi_in_range(v->phi-ve_list[i+1].phi));
276 
277  // two sources of error are (GPS CCN29-19) epsilon/(R sin theta)
278  // and sqrt(2*epsilon/(R (1-cos theta))) and the way things work
279  // out, it is the _smaller_ of the two that is relevant [NB have
280  // changed definition of theta here relative to that used in
281  // CCN29] [NB2: write things so as to avoid zero denominators and
282  // to minimize the multiplications, divisions and above all sqrts
283  // -- that means that c & s are defined including a factor of VR2]
284  c = dot_product(OP,OC);
285  s = fabs(cross_product(OP,OC));
286  double inv_err1 = s * inv_R_EPS_COCIRC;
287  double inv_err2_sq = (R2-c) * inv_R_2EPS_COCIRC;
288  ve_list[i].cocircular_range = pow2(inv_err1) > inv_err2_sq ?
289  1.0/inv_err1 :
290  sqrt(1.0/inv_err2_sq);
291  ve_list[i+1].cocircular_range = ve_list[i].cocircular_range;
292  }
293 }
294 
295 }
The SISCone project has been developed by Gavin Salam and Gregory Soyez
Documentation generated on Mon May 6 2013 11:30:35 for SISCone by  Doxygen 1.8.3.1