/////////////////////////////////////////////////////////////////////////////// // File: quadtree.cpp // // Description: source file for quadtree management (Cquadtree class) // // This file is part of the SISCone project. // // For more details, see http://projects.hepforge.org/siscone // // // // Copyright (c) 2006 Gavin Salam and Gregory Soyez // // // // This program is free software; you can redistribute it and/or modify // // it under the terms of the GNU General Public License as published by // // the Free Software Foundation; either version 2 of the License, or // // (at your option) any later version. // // // // This program is distributed in the hope that it will be useful, // // but WITHOUT ANY WARRANTY; without even the implied warranty of // // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the // // GNU General Public License for more details. // // // // You should have received a copy of the GNU General Public License // // along with this program; if not, write to the Free Software // // Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301, USA // // // // $Revision:: 859 $// // $Date:: 2012-11-28 02:49:23 +0100 (Wed, 28 Nov 2012) $// /////////////////////////////////////////////////////////////////////////////// #include "quadtree.h" #include #include #include namespace siscone{ using namespace std; /******************************************************************* * Cquadtree implementation * * Implementation of a 2D quadtree. * * This class implements the traditional two-dimensional quadtree. * * The elements at each node are of 'Cmomentum' type. * *******************************************************************/ // default ctor //-------------- Cquadtree::Cquadtree(){ v = NULL; children[0][0] = children[0][1] = children[1][0] = children[1][1] = NULL; has_child = false; } // ctor with initialisation (see init for details) //-------------------------- Cquadtree::Cquadtree(double _x, double _y, double _half_size_x, double _half_size_y){ v = NULL; children[0][0] = children[0][1] = children[1][0] = children[1][1] = NULL; has_child = false; init(_x, _y, _half_size_x, _half_size_y); } // default destructor // at destruction, everything is destroyed except // physical values at the leaves //------------------------------------------------ Cquadtree::~Cquadtree(){ if (has_child){ if (v!=NULL) delete v; delete children[0][0]; delete children[0][1]; delete children[1][0]; delete children[1][1]; } } /* * init the tree. * By initializing the tree, we mean setting the cell parameters * and preparing the object to act as a seed for a new tree. * - _x x-position of the center * - _y y-position of the center * - half_size_x half x-size of the cell * - half_size_y half y-size of the cell * return 0 on success, 1 on error. Note that if the cell * is already filled, we return an error. ******************************************************************/ int Cquadtree::init(double _x, double _y, double _half_size_x, double _half_size_y){ if (v!=NULL) return 1; centre_x = _x; centre_y = _y; half_size_x = _half_size_x; half_size_y = _half_size_y; return 0; } /* * adding a particle to the tree. * This method adds one vector to the quadtree structure which * is updated consequently. * - v vector to add * return 0 on success 1 on error ******************************************************************/ int Cquadtree::add(Cmomentum *v_add){ // Description of the method: // -------------------------- // the addition process goes as follows: // 1. check if the cell is empty, in which case, add the particle // here and leave. // 2. If there is a unique particle already inside, // (a) create children // (b) forward the existing particle to the appropriate child // 3. Add current particle to this cell and forward to the // adequate child // NOTE: we assume in the whole procedure that the particle is // indeed inside the cell ! // step 1: the case of empty cells if (v==NULL){ v = v_add; return 0; } // step 2: additional work if 1! particle already present // we use the fact that only 1-particle systems have no child if (!has_child){ double new_half_size_x = 0.5*half_size_x; double new_half_size_y = 0.5*half_size_y; // create children children[0][0] = new Cquadtree(centre_x-new_half_size_x, centre_y-new_half_size_y, new_half_size_x, new_half_size_y); children[0][1] = new Cquadtree(centre_x-new_half_size_x, centre_y+new_half_size_y, new_half_size_x, new_half_size_y); children[1][0] = new Cquadtree(centre_x+new_half_size_x, centre_y-new_half_size_y, new_half_size_x, new_half_size_y); children[1][1] = new Cquadtree(centre_x+new_half_size_x, centre_y+new_half_size_y, new_half_size_x, new_half_size_y); has_child = true; // forward to child //? The following line assumes 'true'==1 and 'false'==0 // Note: v being a single particle, eta and phi are correct children[v->eta>centre_x][v->phi>centre_y]->add(v); // copy physical params v = new Cmomentum(*v); } // step 3: add new particle // Note: v_add being a single particle, eta and phi are correct children[v_add->eta>centre_x][v_add->phi>centre_y]->add(v_add); *v+=*v_add; return 0; } /* * circle intersection. * computes the intersection with a circle of given centre and radius. * The output takes the form of a quadtree with all squares included * in the circle. * - cx circle centre x coordinate * - cy circle centre y coordinate * - cR2 circle radius SQUARED * return the checksum for the intersection ******************************************************************/ Creference Cquadtree::circle_intersect(double cx, double cy, double cR2){ // Description of the method: // -------------------------- // 1. check if cell is empty => no intersection // 2. if cell has 1! particle, check if it is inside the circle. // If yes, add it and return, if not simply return. // 3. check if the circle intersects the square. If not, return. // 4. check if the square is inside the circle. // If yes, add it to qt and return. // 5. check intersections with children. // step 1: if there is no particle inside te square, no reason to go further if (v==NULL) return Creference(); double dx, dy; // step 2: if there is only one particle inside the square, test if it is in // the circle, in which case return associated reference if (!has_child){ // compute the distance // Note: v has only one particle => eta and phi are defined dx = cx - v->eta; dy = fabs(cy - v->phi); if (dy>M_PI) dy -= 2.0*M_PI; // test distance if (dx*dx+dy*dyref; } return Creference(); } // step 3: check if there is an intersection //double ryp, rym; double dx_c, dy_c; // store distance with the centre of the square dx_c = fabs(cx-centre_x); dy_c = fabs(cy-centre_y); if (dy_c>M_PI) dy_c = 2.0*M_PI-dy_c; // compute (minimal) the distance (pay attention to the periodicity in phi). dx = dx_c-half_size_x; if (dx<0) dx=0; dy = dy_c-half_size_y; if (dy<0) dy=0; // check the distance if (dx*dx+dy*dy>=cR2){ // no intersection return Creference(); } // step 4: check if included // compute the (maximal) distance dx = dx_c+half_size_x; dy = dy_c+half_size_y; if (dy>M_PI) dy = M_PI; // compute the distance if (dx*dx+dy*dyref; } // step 5: the square is not fully in. Recurse to children return children[0][0]->circle_intersect(cx, cy, cR2) + children[0][1]->circle_intersect(cx, cy, cR2) + children[1][0]->circle_intersect(cx, cy, cR2) + children[1][1]->circle_intersect(cx, cy, cR2); } /* * output a data file for drawing the grid. * This can be used to output a data file containing all the * grid subdivisions. The file contents is as follows: * first and second columns give center of the cell, the third * gives the size. * - flux opened stream to write to * return 0 on success, 1 on error ******************************************************************/ int Cquadtree::save(FILE *flux){ if (flux==NULL) return 1; if (has_child){ fprintf(flux, "%e\t%e\t%e\t%e\n", centre_x, centre_y, half_size_x, half_size_y); children[0][0]->save(flux); children[0][1]->save(flux); children[1][0]->save(flux); children[1][1]->save(flux); } return 0; } /* * output a data file for drawing the tree leaves. * This can be used to output a data file containing all the * tree leaves. The file contents is as follows: * first and second columns give center of the cell, the third * gives the size. * - flux opened stream to write to * return 0 on success, 1 on error ******************************************************************/ int Cquadtree::save_leaves(FILE *flux){ if (flux==NULL) return 1; if (has_child){ if (children[0][0]!=NULL) children[0][0]->save_leaves(flux); if (children[0][1]!=NULL) children[0][1]->save_leaves(flux); if (children[1][0]!=NULL) children[1][0]->save_leaves(flux); if (children[1][1]!=NULL) children[1][1]->save_leaves(flux); } else { fprintf(flux, "%e\t%e\t%e\t%e\n", centre_x, centre_y, half_size_x, half_size_y); } return 0; } }