SISCone  3.0.5
geom_2d.cpp
1 // File: geom_2d.cpp //
3 // Description: source file for two-dimensional geometry tools //
4 // This file is part of the SISCone project. //
5 // WARNING: this is not the main SISCone trunk but //
6 // an adaptation to spherical coordinates //
7 // For more details, see http://projects.hepforge.org/siscone //
8 // //
9 // Copyright (c) 2006-2008 Gavin Salam and Gregory Soyez //
10 // //
11 // This program is free software; you can redistribute it and/or modify //
12 // it under the terms of the GNU General Public License as published by //
13 // the Free Software Foundation; either version 2 of the License, or //
14 // (at your option) any later version. //
15 // //
16 // This program is distributed in the hope that it will be useful, //
17 // but WITHOUT ANY WARRANTY; without even the implied warranty of //
18 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the //
19 // GNU General Public License for more details. //
20 // //
21 // You should have received a copy of the GNU General Public License //
22 // along with this program; if not, write to the Free Software //
23 // Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301, USA //
24 // //
25 // $Revision:: 378 $//
26 // $Date:: 2016-02-24 15:10:38 +0100 (Wed, 24 Feb 2016) $//
28 
29 #include "geom_2d.h"
30 #include <algorithm>
31 
32 namespace siscone_spherical{
33 
34 #define PHI_RANGE_MASK 0xFFFFFFFF
35 
36 /*********************************************************
37  * class CSphtheta_phi_range implementation *
38  * class for holding a covering range in eta-phi *
39  * *
40  * This class deals with ranges in the eta-phi plane. It *
41  * implements methods to test if two ranges overlap and *
42  * to take the union of two overlapping intervals. *
43  *********************************************************/
44 
45 using namespace std;
46 
47 // static member default init
48 //----------------------------
50 double CSphtheta_phi_range::theta_max = M_PI;
51 
52 // default ctor
53 //--------------
55  theta_range = 0;
56  phi_range = 0;
57 }
58 
59 // ctor with initialisation
60 // we initialise with a centre (in eta,phi) and a radius
61 // - c_theta theta coordinate of the centre
62 // - c_phi phi coordinate of the centre
63 // - R radius
64 //-------------------------------------------------------
65 CSphtheta_phi_range::CSphtheta_phi_range(double c_theta, double c_phi, double R){
66  // determination of the eta range
67  //-------------------------------
68  double xmin = max(c_theta-R,theta_min+0.00001);
69  double xmax = min(c_theta+R,theta_max-0.00001);
70 
71  unsigned int cell_min = get_theta_cell(xmin);
72  unsigned int cell_max = get_theta_cell(xmax);
73 
74  // warning: if cell_max==2^31, 2*cell_max==0 hence,
75  // even if the next formula is formally (2*cell_max-cell_min),
76  // expressing it as (cell_max-cell_min)+cell_max is safe.
77  theta_range = (cell_max-cell_min)+cell_max;
78 
79  // determination of the phi range
80  // !! taking care of periodicity !!
81  // !! and the theta dependence !!
82  //---------------------------------
83  double ymin,ymax;
84  double extra = asin(R/M_PI);
85  // if the theta range comes too close to the endpoints (theta=0 or
86  // theta=pi), then keep the full phi range
87  if (xmin<=theta_min+extra){
88  ymin = -M_PI+0.00001;
89  ymax = M_PI-0.00001;
90  } else if (xmax>=theta_max-extra){
91  ymin = -M_PI+0.00001;
92  ymax = M_PI-0.00001;
93  } else {
94  extra = max(1.0/sin(xmin), 1.0/sin(xmax));
95  ymin = (c_phi-R)*extra;
96  while (ymin<-M_PI) ymin+=twopi;
97  while (ymin> M_PI) ymin-=twopi;
98  ymax = (c_phi-R)*extra;
99  while (ymax<-M_PI) ymax+=twopi;
100  while (ymax> M_PI) ymax-=twopi;
101  }
102  cell_min = get_phi_cell(ymin);
103  cell_max = get_phi_cell(ymax);
104 
105  // Also, if the interval goes through pi, inversion is needed
106  if (ymax>ymin)
107  phi_range = (cell_max-cell_min)+cell_max;
108  else {
109  phi_range = (cell_min==cell_max)
110  ? PHI_RANGE_MASK
111  : ((PHI_RANGE_MASK^(cell_min-cell_max)) + cell_max);
112  }
113 }
114 
115 // assignment of range
116 // - r range to assign to current one
117 //---------------------------------------
119  theta_range = r.theta_range;
120  phi_range = r.phi_range;
121 
122  return *this;
123 }
124 
125 // add a particle to the range
126 // - eta eta coordinate of the particle
127 // - phi phi coordinate of the particle
128 // \return 0 on success, 1 on error
129 //----------------------------------------
130 int CSphtheta_phi_range::add_particle(const double theta, const double phi){
131  // get the theta cell
132  unsigned int theta_cell = get_theta_cell(theta);
133 
134  // deal with the eta coordinate
135  theta_range |= theta_cell;
136 
137  // deal with the phi coordinate
138  //
139  // watch out: if the theta_cell includes theta==0 or theta==pi,
140  // incude the full phi range
141  if ((theta_cell == 0x1) || (theta_cell == 0x80000000))
142  phi_range = 0xffffffff;
143  else
144  phi_range |= get_phi_cell(phi);
145 
146  return 0;
147 }
148 
149 
150 // test overlap
151 // - r1 first range
152 // - r2 second range
153 // return true if overlap, false otherwise.
154 //------------------------------------------
155 bool is_range_overlap(const CSphtheta_phi_range &r1, const CSphtheta_phi_range &r2){
156  // check overlap in eta AND phi
157  return ((r1.theta_range & r2.theta_range) && (r1.phi_range & r2.phi_range));
158 }
159 
160 // compute union
161 // Note: we assume that the two intervals overlap
162 // - r1 first range
163 // - r2 second range
164 // \return union of the two ranges
165 //------------------------------------------
166 const CSphtheta_phi_range range_union (const CSphtheta_phi_range &r1, const CSphtheta_phi_range &r2){
167  CSphtheta_phi_range tmp;
168 
169  // compute union in eta
170  tmp.theta_range = r1.theta_range | r2.theta_range;
171 
172  // compute union in phi
173  tmp.phi_range = r1.phi_range | r2.phi_range;
174 
175  return tmp;
176 }
177 
178 }
unsigned int phi_range
phi range as a binary coding of covered cells
Definition: geom_2d.h:146
class for holding a covering range in eta-phi
Definition: geom_2d.h:52
int add_particle(const double theta, const double phi)
add a particle to the range
Definition: geom_2d.cpp:130
CSphtheta_phi_range & operator=(const CSphtheta_phi_range &r)
assignment of range
Definition: geom_2d.cpp:118
static double theta_max
maximal value for theta (set to pi)
Definition: geom_2d.h:82
unsigned int theta_range
theta range as a binary coding of covered cells
Definition: geom_2d.h:75
static double theta_min
extremal value for theta
Definition: geom_2d.h:81
unsigned int phi_range
phi range as a binary coding of covered cells
Definition: geom_2d.h:78
const double twopi
definition of 2*M_PI which is useful a bit everyhere!
Definition: defines.h:114