source: HiSusy/trunk/Delphes/Delphes-3.0.9/external/fastjet/plugins/SISCone/vicinity.cc @ 5

Last change on this file since 5 was 5, checked in by zerwas, 11 years ago

update to Delphes-3.0.9

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 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
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.