Fork me on GitHub

source: git/external/fastjet/plugins/SISCone/vicinity.cc

Last change on this file was 1d208a2, checked in by Pavel Demin <pavel.demin@…>, 8 years ago

update FastJet library to 3.2.1 and Nsubjettiness library to 2.2.4

  • Property mode set to 100644
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:: 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
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 ve_list = NULL;
83 set_particle_list(_particle_list);
84}
85
86// default destructor
87//--------------------
88Cvicinity::~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 ************************************************************/
103void 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 ************************************************************/
176void 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
208inline 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 **********************************************************/
220void 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}
Note: See TracBrowser for help on using the repository browser.