source: HiSusy/trunk/Delphes/Delphes-3.0.9/external/fastjet/plugins/SISCone/quadtree.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: 10.1 KB
Line 
1///////////////////////////////////////////////////////////////////////////////
2// File: quadtree.cpp                                                        //
3// Description: source file for quadtree management (Cquadtree 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 "quadtree.h"
28#include <math.h>
29#include <stdio.h>
30#include <iostream>
31
32namespace siscone{
33
34using namespace std;
35
36/*******************************************************************
37 * Cquadtree implementation                                        *
38 * Implementation of a 2D quadtree.                                *
39 * This class implements the traditional two-dimensional quadtree. *
40 * The elements at each node are of 'Cmomentum' type.              *
41 *******************************************************************/
42
43// default ctor
44//--------------
45Cquadtree::Cquadtree(){
46  v = NULL;
47
48  children[0][0] = children[0][1] = children[1][0] = children[1][1] = NULL;
49  has_child = false;
50}
51
52
53// ctor with initialisation (see init for details)
54//--------------------------
55Cquadtree::Cquadtree(double _x, double _y, double _half_size_x, double _half_size_y){
56  v = NULL;
57
58  children[0][0] = children[0][1] = children[1][0] = children[1][1] = NULL;
59  has_child = false;
60
61  init(_x, _y, _half_size_x, _half_size_y);
62}
63
64
65// default destructor
66// at destruction, everything is destroyed except
67// physical values at the leaves
68//------------------------------------------------
69Cquadtree::~Cquadtree(){
70  if (has_child){
71    if (v!=NULL) delete v;
72    delete children[0][0];
73    delete children[0][1];
74    delete children[1][0];
75    delete children[1][1];
76  }
77}
78
79
80/*
81 * init the tree.
82 * By initializing the tree, we mean setting the cell parameters
83 * and preparing the object to act as a seed for a new tree.
84 *  - _x           x-position of the center
85 *  - _y           y-position of the center
86 *  - half_size_x  half x-size of the cell
87 *  - half_size_y  half y-size of the cell
88 * return 0 on success, 1 on error. Note that if the cell
89 *        is already filled, we return an error.
90 ******************************************************************/
91int Cquadtree::init(double _x, double _y, double _half_size_x, double _half_size_y){
92  if (v!=NULL)
93    return 1;
94
95  centre_x = _x;
96  centre_y = _y;
97  half_size_x = _half_size_x;
98  half_size_y = _half_size_y;
99
100  return 0;
101}
102
103
104/*
105 * adding a particle to the tree.
106 * This method adds one vector to the quadtree structure which
107 * is updated consequently.
108 *  - v   vector to add
109 * return 0 on success 1 on error
110 ******************************************************************/
111int Cquadtree::add(Cmomentum *v_add){
112  // Description of the method:
113  // --------------------------
114  // the addition process goes as follows:
115  //  1. check if the cell is empty, in which case, add the particle
116  //     here and leave.
117  //  2. If there is a unique particle already inside,
118  //      (a) create children
119  //      (b) forward the existing particle to the appropriate child
120  //  3. Add current particle to this cell and forward to the
121  //     adequate child
122  // NOTE: we assume in the whole procedure that the particle is
123  //       indeed inside the cell !
124
125  // step 1: the case of empty cells
126  if (v==NULL){
127    v = v_add;
128    return 0;
129  }
130
131  // step 2: additional work if 1! particle already present
132  //         we use the fact that only 1-particle systems have no child
133  if (!has_child){
134    double new_half_size_x = 0.5*half_size_x;
135    double new_half_size_y = 0.5*half_size_y;
136    // create children
137    children[0][0] = new Cquadtree(centre_x-new_half_size_x, centre_y-new_half_size_y,
138                                   new_half_size_x, new_half_size_y);
139    children[0][1] = new Cquadtree(centre_x-new_half_size_x, centre_y+new_half_size_y,
140                                   new_half_size_x, new_half_size_y);
141    children[1][0] = new Cquadtree(centre_x+new_half_size_x, centre_y-new_half_size_y,
142                                   new_half_size_x, new_half_size_y);
143    children[1][1] = new Cquadtree(centre_x+new_half_size_x, centre_y+new_half_size_y,
144                                   new_half_size_x, new_half_size_y);
145
146    has_child = true;
147
148    // forward to child
149    //? The following line assumes 'true'==1 and 'false'==0
150    // Note: v being a single particle, eta and phi are correct
151    children[v->eta>centre_x][v->phi>centre_y]->add(v);
152
153    // copy physical params
154    v = new Cmomentum(*v);
155  }
156
157  // step 3: add new particle
158  // Note: v_add being a single particle, eta and phi are correct
159  children[v_add->eta>centre_x][v_add->phi>centre_y]->add(v_add);
160  *v+=*v_add;
161
162  return 0;
163}
164
165
166/*
167 * circle intersection.
168 * computes the intersection with a circle of given centre and radius.
169 * The output takes the form of a quadtree with all squares included
170 * in the circle.
171 *  - cx    circle centre x coordinate
172 *  - cy    circle centre y coordinate
173 *  - cR2   circle radius SQUARED
174 * return the checksum for the intersection
175 ******************************************************************/
176Creference Cquadtree::circle_intersect(double cx, double cy, double cR2){
177  // Description of the method:
178  // --------------------------
179  // 1. check if cell is empty => no intersection
180  // 2. if cell has 1! particle, check if it is inside the circle.
181  //    If yes, add it and return, if not simply return.
182  // 3. check if the circle intersects the square. If not, return.
183  // 4. check if the square is inside the circle.
184  //    If yes, add it to qt and return.
185  // 5. check intersections with children.
186
187  // step 1: if there is no particle inside te square, no reason to go further
188  if (v==NULL)
189    return Creference();
190
191  double dx, dy;
192
193  // step 2: if there is only one particle inside the square, test if it is in
194  //         the circle, in which case return associated reference
195  if (!has_child){
196    // compute the distance
197    // Note: v has only one particle => eta and phi are defined
198    dx = cx - v->eta;
199    dy = fabs(cy - v->phi);
200    if (dy>M_PI) 
201      dy -= 2.0*M_PI;
202
203    // test distance
204    if (dx*dx+dy*dy<cR2){
205      return v->ref;
206    }
207
208    return Creference();
209  }
210
211  // step 3: check if there is an intersection
212  //double ryp, rym;
213  double dx_c, dy_c;
214
215  // store distance with the centre of the square
216  dx_c = fabs(cx-centre_x);
217  dy_c = fabs(cy-centre_y);
218  if (dy_c>M_PI) dy_c = 2.0*M_PI-dy_c;
219
220  // compute (minimal) the distance (pay attention to the periodicity in phi).
221  dx = dx_c-half_size_x;
222  if (dx<0) dx=0;
223  dy = dy_c-half_size_y;
224  if (dy<0) dy=0;
225
226  // check the distance
227  if (dx*dx+dy*dy>=cR2){
228    // no intersection
229    return Creference();
230  }
231
232  // step 4: check if included
233
234  // compute the (maximal) distance
235  dx = dx_c+half_size_x;
236  dy = dy_c+half_size_y;
237  if (dy>M_PI) dy = M_PI;
238
239  // compute the distance
240  if (dx*dx+dy*dy<cR2){
241    return v->ref;
242  }
243
244  // step 5: the square is not fully in. Recurse to children
245  return children[0][0]->circle_intersect(cx, cy, cR2)
246    + children[0][1]->circle_intersect(cx, cy, cR2)
247    + children[1][0]->circle_intersect(cx, cy, cR2)
248    + children[1][1]->circle_intersect(cx, cy, cR2);
249}
250
251
252/*
253 * output a data file for drawing the grid.
254 * This can be used to output a data file containing all the
255 * grid subdivisions. The file contents is as follows:
256 * first and second columns give center of the cell, the third
257 * gives the size.
258 *  - flux  opened stream to write to
259 * return 0 on success, 1 on error
260 ******************************************************************/
261int Cquadtree::save(FILE *flux){
262
263  if (flux==NULL)
264    return 1;
265
266  if (has_child){
267    fprintf(flux, "%e\t%e\t%e\t%e\n", centre_x, centre_y, half_size_x, half_size_y);
268    children[0][0]->save(flux);
269    children[0][1]->save(flux);
270    children[1][0]->save(flux);
271    children[1][1]->save(flux);
272  }
273
274  return 0;
275}
276
277
278/*
279 * output a data file for drawing the tree leaves.
280 * This can be used to output a data file containing all the
281 * tree leaves. The file contents is as follows:
282 * first and second columns give center of the cell, the third
283 * gives the size.
284 *  - flux  opened stream to write to
285 * return 0 on success, 1 on error
286 ******************************************************************/
287int Cquadtree::save_leaves(FILE *flux){
288
289  if (flux==NULL)
290    return 1;
291
292  if (has_child){
293    if (children[0][0]!=NULL) children[0][0]->save_leaves(flux);
294    if (children[0][1]!=NULL) children[0][1]->save_leaves(flux);
295    if (children[1][0]!=NULL) children[1][0]->save_leaves(flux);
296    if (children[1][1]!=NULL) children[1][1]->save_leaves(flux);
297  } else {
298    fprintf(flux, "%e\t%e\t%e\t%e\n", centre_x, centre_y, half_size_x, half_size_y);
299  }
300
301  return 0;
302}
303
304}
Note: See TracBrowser for help on using the repository browser.