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 // WARNING: this is not the main SISCone trunk but //
6 // an adaptation to spherical coordinates //
7 // For more details, see http://projects.hepforge.org/siscone //
8 // //
9 // Copyright (c) 2006-2008 Gavin Salam and Gregory Soyez //
10 // //
11 // This program is free software; you can redistribute it and/or modify //
12 // it under the terms of the GNU General Public License as published by //
13 // the Free Software Foundation; either version 2 of the License, or //
14 // (at your option) any later version. //
15 // //
16 // This program is distributed in the hope that it will be useful, //
17 // but WITHOUT ANY WARRANTY; without even the implied warranty of //
18 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the //
19 // GNU General Public License for more details. //
20 // //
21 // You should have received a copy of the GNU General Public License //
22 // along with this program; if not, write to the Free Software //
23 // Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301, USA //
24 // //
25 // $Revision:: 255 $//
26 // $Date:: 2008-07-12 17:40:35 +0200 (Sat, 12 Jul 2008) $//
28 
29 #include "vicinity.h"
30 #include <math.h>
31 #include <algorithm>
32 #include <iostream>
33 
34 namespace siscone_spherical{
35 
36 using namespace std;
37 
38 /*************************************************************
39  * CSphvicinity_elm implementation *
40  * element in the vicinity of a parent. *
41  * class used to manage one points in the vicinity *
42  * of a parent point. *
43  *************************************************************/
44 
45 // ordering pointers to CSphvicinity_elm
46 //---------------------------------------
47 bool ve_less(CSphvicinity_elm *ve1, CSphvicinity_elm *ve2){
48  return ve1->angle < ve2->angle;
49 }
50 
51 
52 /*************************************************************
53  * CSphvicinity implementation *
54  * list of element in the vicinity of a parent. *
55  * class used to manage the points which are in the vicinity *
56  * of a parent point. The construction of the list can be *
57  * made from a list of points or from a quadtree. *
58  *************************************************************/
59 
60 // default constructor
61 //---------------------
63  n_part = 0;
64 
65  ve_list = NULL;
66 #ifdef USE_QUADTREE_FOR_STABILITY_TEST
67  quadtree = NULL;
68 #endif
69 
70  parent = NULL;
71  VR2 = VR = 0.0;
72 
73 }
74 
75 // constructor with initialisation
76 //---------------------------------
77 CSphvicinity::CSphvicinity(vector<CSphmomentum> &_particle_list){
78  parent = NULL;
79  ve_list = NULL;
80 #ifdef USE_QUADTREE_FOR_STABILITY_TEST
81  quadtree = NULL;
82 #endif
83  cosVR = VR2 = tan2R = VR = 0.0;
84 
85  set_particle_list(_particle_list);
86 }
87 
88 // default destructor
89 //--------------------
91  if (ve_list!=NULL)
92  delete[] ve_list;
93 
94 #ifdef USE_QUADTREE_FOR_STABILITY_TEST
95  if (quadtree!=NULL)
96  delete quadtree;
97 #endif
98 }
99 
100 /*
101  * set the particle_list
102  * - particle_list list of particles (type CSphmomentum)
103  * - n number of particles in the list
104  ************************************************************/
105 void CSphvicinity::set_particle_list(vector<CSphmomentum> &_particle_list){
106  int i,j;
107 #ifdef USE_QUADTREE_FOR_STABILITY_TEST
108  double eta_max=0.0;
109 #endif
110 
111  // if the particle list is not empty, destroy it !
112  if (ve_list!=NULL){
113  delete[] ve_list;
114  }
115  vicinity.clear();
116 #ifdef USE_QUADTREE_FOR_STABILITY_TEST
117  if (quadtree!=NULL)
118  delete quadtree;
119 #endif
120 
121  // allocate memory array for particles
122  // Note: - we compute max for |eta|
123  // - we allocate indices to particles
124  n_part = 0;
125  plist.clear();
126  pincluded.clear();
127  for (i=0;i<(int) _particle_list.size();i++){
128  // if a particle is colinear with the beam (infinite rapidity)
129  // we do not take it into account
130  //if (fabs(_particle_list[i].pz)!=_particle_list[i].E){
131  plist.push_back(_particle_list[i]);
132  pincluded.push_back(siscone::Cvicinity_inclusion()); // zero inclusion status
133 
134  // the parent_index is handled in the split_merge because
135  // of our multiple-pass procedure.
136  // Hence, it is not required here any longer.
137  // plist[n_part].parent_index = i;
138  plist[n_part].index = n_part;
139 
140  // make sure the reference is randomly created
141  plist[n_part].ref.randomize();
142 
143 #ifdef USE_QUADTREE_FOR_STABILITY_TEST
144  if (fabs(plist[n_part].eta)>eta_max) eta_max=fabs(plist[n_part].eta);
145 #endif
146  n_part++;
147  //}
148  }
149 
150  // allocate quadtree and vicinity_elm list
151  // note: we set phi in [-pi:pi] as it is the natural range for atan2!
152  ve_list = new CSphvicinity_elm[2*n_part];
153 #ifdef USE_QUADTREE_FOR_STABILITY_TEST
154  eta_max+=0.1;
155  quadtree = new siscone::Cquadtree(0.0, 0.0, eta_max, M_PI);
156 #endif
157 
158  // append particle to the vicinity_elm list
159  j = 0;
160  for (i=0;i<n_part;i++){
161 #ifdef USE_QUADTREE_FOR_STABILITY_TEST
162  quadtree->add(&plist[i]);
163 #endif
164  ve_list[j].v = ve_list[j+1].v = &plist[i];
165  ve_list[j].is_inside = ve_list[j+1].is_inside = &(pincluded[i]);
166  j+=2;
167  }
168 
169 }
170 
171 
172 /*
173  * build the vicinity list from a list of points.
174  * - _parent reference particle
175  * - _VR vicinity radius
176  ************************************************************/
177 void CSphvicinity::build(CSphmomentum *_parent, double _VR){
178  int i;
179 
180  // set parent and radius
181  parent = _parent;
182 
183  VR = _VR;
184  VR2 = VR*VR;
185  cosVR = cos(VR);
186  R2 = 0.25*VR2;
187  R = 0.5*VR;
188  double tmp = tan(R);
189  tan2R = tmp*tmp;
190 
191  D2_R = 2.0*(1-cos(R));
192  //tmp = sqrt(D2_R);
193  inv_R_EPS_COCIRC = 1.0 / R / EPSILON_COCIRCULAR;
194  inv_R_2EPS_COCIRC = 0.5 / R / EPSILON_COCIRCULAR;
195 
196  // clear vicinity
197  vicinity.clear();
198 
199  // init parent variables
200  // we cpte the direction of the centre and two orthogonal ones
201  // to measure the angles. Those are taken orthogonal to the
202  // axis of smallest components (of the centre) to increase precision
203  parent_centre = (*parent)/parent->_norm;
204  parent_centre.get_angular_directions(angular_dir1, angular_dir2);
205  angular_dir1 /= angular_dir1._norm;
206  angular_dir2 /= angular_dir2._norm;
207 
208  // really browse the particle list
209  for (i=0;i<n_part;i++){
210  append_to_vicinity(&plist[i]);
211  }
212 
213  // sort the vicinity
214  sort(vicinity.begin(), vicinity.end(), ve_less);
215 
216  vicinity_size = vicinity.size();
217 }
218 
219 
221 //TODO//
222 inline double sort_angle(double s, double c){
223  if (s==0) return (c>0) ? 0.0 : 2.0;
224  double t=c/s;
225  return (s>0) ? 1-t/(1+fabs(t)) : 3-t/(1+fabs(t));
226 }
227 
228 
229 /*
230  * append a particle to the 'vicinity' list after
231  * having computed the angular-ordering quantities
232  * - v vector to test
233  **********************************************************/
235  // skip the particle itself)
236  if (v==parent)
237  return;
238 
239  int i=2*(v->index);
240 
241  // compute the distance of the i-th particle with the parent
242  double dot = dot_product3(parent_centre,*v);
243  CSph3vector vnormal = *v;
244  vnormal/=v->_norm;
245  dot/=v->_norm;
246 
247  // really check if the distance is less than VR
248  if (dot>cosVR){
249  CSph3vector cross = cross_product3(parent_centre,vnormal);
250 
251  // for the centres
252  CSph3vector median = (parent_centre+vnormal);
253  double amplT = sqrt((tan2R*(1+dot)+(dot-1))*(1+dot));
254  CSph3vector transverse = amplT*cross/cross._norm;
255 
256  // first angle (+)
257  ve_list[i].centre = median + transverse;
258  ve_list[i].centre.build_norm();
259  ve_list[i].centre/=ve_list[i].centre._norm;
260  CSph3vector diff = ve_list[i].centre - parent_centre;
261  //ve_list[i].angle = atan2(dot_product3(angular_dir2, diff),dot_product3(angular_dir1, diff));
262  ve_list[i].angle = sort_angle(dot_product3(angular_dir2, diff),dot_product3(angular_dir1, diff));
263  ve_list[i].side = true;
264  ve_list[i].cocircular.clear();
265  vicinity.push_back(&(ve_list[i]));
266 
267  // second angle (-)
268  ve_list[i+1].centre = median - transverse;
269  ve_list[i+1].centre.build_norm();
270  ve_list[i+1].centre/=ve_list[i+1].centre._norm;
271  diff = ve_list[i+1].centre - parent_centre;
272  ve_list[i+1].angle = sort_angle(dot_product3(angular_dir2, diff),dot_product3(angular_dir1, diff));
273  ve_list[i+1].side = false;
274  ve_list[i+1].cocircular.clear();
275  vicinity.push_back(&(ve_list[i+1]));
276 
277  // now work out the cocircularity range for the two points (range
278  // of angle within which the points stay within a distance
279  // EPSILON_COCIRCULAR of circule
280  // P = parent; C = child; O = Origin (center of circle)
281  CSph3vector OP = parent_centre - ve_list[i+1].centre;
282  CSph3vector OC = vnormal - ve_list[i+1].centre;
283 
284  // two sources of error are (GPS CCN29-19) epsilon/(R sin theta)
285  // and sqrt(2*epsilon/(R (1-cos theta))) and the way things work
286  // out, it is the _smaller_ of the two that is relevant [NB have
287  // changed definition of theta here relative to that used in
288  // CCN29] [NB2: write things so as to avoid zero denominators and
289  // to minimize the multiplications, divisions and above all sqrts
290  // -- that means that c & s are defined including a factor of VR2]
291  double inv_err1 = cross_product3(OP,OC)._norm * inv_R_EPS_COCIRC;
292  double inv_err2_sq = (D2_R-dot_product3(OP,OC)) * inv_R_2EPS_COCIRC;
293  ve_list[i].cocircular_range = siscone::pow2(inv_err1) > inv_err2_sq ?
294  1.0/inv_err1 :
295  sqrt(1.0/inv_err2_sq);
296  ve_list[i+1].cocircular_range = ve_list[i].cocircular_range;
297  }
298 }
299 
300 }
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