Fork me on GitHub

source: svn/trunk/external/fastjet/plugins/SISCone/vicinity.cc@ 1397

Last change on this file since 1397 was 859, checked in by Pavel Demin, 12 years ago

update fastjet to version 3.0.3

  • Property svn:eol-style set to native
  • Property svn:keywords set to Id Revision Date
File size: 9.7 KB
Line 
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 01:49:23 +0000 (Wed, 28 Nov 2012) $//
25///////////////////////////////////////////////////////////////////////////////
26
27#include "vicinity.h"
28#include <math.h>
29#include <algorithm>
30#include <iostream>
31
32namespace siscone{
33
34using 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//------------------------------------
45bool 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//---------------------
60Cvicinity::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//---------------------------------
75Cvicinity::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//--------------------
87Cvicinity::~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 ************************************************************/
102void 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 ************************************************************/
175void 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
207inline 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 **********************************************************/
219void 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}
Note: See TracBrowser for help on using the repository browser.