| 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:: 123 $//
|
|---|
| 24 | // $Date:: 2007-03-01 02:52:16 +0100 (Thu, 01 Mar 2007) $//
|
|---|
| 25 | ///////////////////////////////////////////////////////////////////////////////
|
|---|
| 26 |
|
|---|
| 27 | #include "vicinity.h"
|
|---|
| 28 | #include <math.h>
|
|---|
| 29 | #include <algorithm>
|
|---|
| 30 | #include <iostream>
|
|---|
| 31 |
|
|---|
| 32 | namespace siscone{
|
|---|
| 33 |
|
|---|
| 34 | using 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 | //------------------------------------
|
|---|
| 45 | bool 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 | //---------------------
|
|---|
| 60 | Cvicinity::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 | //---------------------------------
|
|---|
| 75 | Cvicinity::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 | //--------------------
|
|---|
| 87 | Cvicinity::~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 | ************************************************************/
|
|---|
| 102 | void 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 | ************************************************************/
|
|---|
| 175 | void 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
|
|---|
| 207 | inline 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 | **********************************************************/
|
|---|
| 219 | void 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 | }
|
|---|
