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:: 388 $//
24 // $Date:: 2016-03-03 10:42:25 +0100 (Thu, 03 Mar 2016) $//
25 ///////////////////////////////////////////////////////////////////////////////
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 //------------------------------------
ve_less(Cvicinity_elm * ve1,Cvicinity_elm * ve2)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 //---------------------
Cvicinity()60 Cvicinity::Cvicinity(){
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 //---------------------------------
Cvicinity(vector<Cmomentum> & _particle_list)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 //--------------------
~Cvicinity()88 Cvicinity::~Cvicinity(){
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 ************************************************************/
set_particle_list(vector<Cmomentum> & _particle_list)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 ************************************************************/
build(Cmomentum * _parent,double _VR)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
207 /// strictly increasing function of the angle
sort_angle(double s,double c)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 **********************************************************/
append_to_vicinity(Cmomentum * v)220 void Cvicinity::append_to_vicinity(Cmomentum *v){
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 }
297