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:: 859 $// |
---|
24 | // $Date:: 2012-11-28 02:49:23 +0100 (Wed, 28 Nov 2012) $// |
---|
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 | //------------------------------------ |
---|
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 | //--------------------- |
---|
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 | //--------------------------------- |
---|
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 | //-------------------- |
---|
87 | Cvicinity::~Cvicinity(){ |
---|
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 | |
---|
206 | /// strictly increasing function of the angle |
---|
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 | **********************************************************/ |
---|
219 | void Cvicinity::append_to_vicinity(Cmomentum *v){ |
---|
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 | } |
---|