ECCE @ EIC Software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
CylindricalHough.h
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file CylindricalHough.h
1 #ifndef __CYLINDRICALHOUGH__
2 #define __CYLINDRICALHOUGH__
3 
4 #include "CircleHough.h"
5 #include <cmath>
6 #include <iostream>
7 
8 
9 #ifndef M_PI
10 #define M_PI 3.14159265358979323846
11 #endif
12 
13 
15 {
16  public:
17  AngleIndexPair(float ang, unsigned int idx) : angle(ang), index(idx)
18  {
19  float twopi = 2.*M_PI;
20  int a = (int)(angle/twopi);
21  angle -= a*twopi;
22  while(angle < 0.){angle += twopi;}
23  while(angle >= twopi){angle -= twopi;}
24  }
26 
27  bool operator<(const AngleIndexPair& other) const
28  {
29  return angle<other.angle;
30  }
31 
32  static float absDiff(float angle1, float angle2)
33  {
34  float diff = ( angle1 - angle2 );
35  while(diff > M_PI){diff -= 2.*M_PI;}
36  while(diff < -M_PI){diff += 2.*M_PI;}
37  diff = fabs(diff);
38  return diff;
39  }
40 
41  float angle;
42  unsigned int index;
43 };
44 
45 
47 {
48  public:
51 
52  void addPair(AngleIndexPair& angind)
53  {
54  sorted=false;
55  vec.push_back(angind);
56  }
57 
58 
59  void getRangeListSimple(float angle, float error, std::vector<AngleIndexPair*>& result)
60  {
61  result.clear();
62 
63  for(unsigned int i=0;i<vec.size();i++)
64  {
65  if(AngleIndexPair::absDiff(angle, vec[i].angle) <= error)
66  {
67  result.push_back(&(vec[i]));
68  }
69  }
70  }
71 
72  void getRangeList(float angle, float error, std::vector<AngleIndexPair*>& result)
73  {
74  float twopi = 2.*M_PI;
75  int a = (int)(angle/twopi);
76  angle -= a*twopi;
77  while(angle < 0.){angle += twopi;}
78  while(angle >= twopi){angle -= twopi;}
79 
80  if(vec.size() <= 4){return getRangeListSimple(angle, error, result);}
81 
82  result.clear();
83 
84  unsigned int closest = findClosest(angle);
85  //first, traverse upward
86  unsigned int current = closest;
87  unsigned int lowest = 0;
88  unsigned int highest = vec.size()-1;
89  while(true)
90  {
91  if(AngleIndexPair::absDiff(angle, vec[current].angle) <= error)
92  {
93  result.push_back(&(vec[current]));
94  current = (current+1)%(vec.size());
95  if(current==closest){break;}
96  }
97  else
98  {
99  break;
100  }
101  }
102 
103  if(closest==0){return;}
104 
105  //now, traverse downward
106  if(current <= closest)
107  {
108  lowest=current;
109  }
110  else
111  {
112  highest=current;
113  }
114  current = closest-1;
115  while(true)
116  {
117  if(AngleIndexPair::absDiff(angle, vec[current].angle) <= error)
118  {
119  result.push_back(&(vec[current]));
120  if( (current==lowest) || (current==highest) ){break;}
121  current = ((current + vec.size()) - 1)%(vec.size());
122  }
123  else
124  {
125  break;
126  }
127  }
128  }
129 
130  private:
131  unsigned int findClosestSimple(float angle, unsigned int lower, unsigned int upper)
132  {
133  unsigned int closest = lower;
134  float diff = AngleIndexPair::absDiff(vec[closest].angle, angle);
135  for(unsigned int i=(lower+1);i<=upper;i++)
136  {
137  float tempdiff = AngleIndexPair::absDiff(vec[i].angle, angle);
138  if( tempdiff < diff )
139  {
140  closest = i;
141  diff = tempdiff;
142  }
143  }
144 
145  return closest;
146  }
147 
148 
149 
150  unsigned int findClosest(float angle)
151  {
152  if(vec.size() <= 4){return findClosestSimple(angle, 0, vec.size()-1);}
153 
154  if(sorted==false)
155  {
156  std::sort(vec.begin(), vec.end());
157  sorted=true;
158  }
159 
160  unsigned int lower = 0;
161  unsigned int upper = vec.size() - 1;
162  unsigned int middle = vec.size()/2;
163  while(true)
164  {
165  if((upper - lower) <= 4){return findClosestSimple(angle, lower, upper);}
166 
167  if(angle <= vec[middle].angle)
168  {
169  upper = middle;
170  middle = (lower + upper)/2;
171  }
172  else
173  {
174  lower = middle;
175  middle = (lower + upper)/2;
176  }
177  }
178  }
179 
180  std::vector<AngleIndexPair> vec;
181  bool sorted;
182 };
183 
184 
185 class CylindricalHough : public CircleHough
186 {
187  public:
188  CylindricalHough(std::vector<float>& detrad, unsigned int inv_radius_nbin, unsigned int center_angle_nbin, unsigned int dca_origin_nbin, CircleResolution& min_resolution, CircleResolution& max_resolution, CircleRange& range, unsigned int z0_nbin, unsigned int theta_nbin, ZResolution& minzres, ZResolution& maxzres, ZRange& zrange, double sxy=70.e-4, double sz=500.e-4);
190 
191  void customFindHelicesInit(std::vector<SimpleHit3D>& hits, unsigned int min_hits, unsigned int max_hits, unsigned int min_zhits, unsigned int max_zhits, double chi2_cut, float xydiffcut, std::vector<SimpleTrack3D>& tracks, unsigned int maxtracks=0);
192 
193  void addHits(unsigned int zlevel, std::vector<SimpleTrack3D>& temptracks, std::vector<SimpleTrack3D>& tracks, std::vector<float>& params, int tracks_per_hit, float z_cut);
194 
195  void init_ZHough(int z0_nbin, unsigned int theta_nbin, ZResolution& minzres, ZResolution& maxzres, ZRange& zrange);
196 
197  void setVertex(double vx, double vy, double vz)
198  {
199  vertex_x = vx;
200  vertex_y = vy;
201  vertex_z = vz;
202  }
203 
204 
205  void setPhiCut(double pc){phicut=pc;}
206 
207  bool intersect_circles(bool hel, double startx, double starty, double rad_det, double rad_trk, double cx, double cy, double& x, double& y);
208 
209  void setLayerResolution(std::vector<double>& lxy, std::vector<double>& lz);
210  void setVertexResolution(double vxy, double vz);
211 
212  private:
213  std::vector<AngleIndexList> angle_list;
214  std::vector<float> detector_radii;
216  double phicut;
217 
218 };
219 
220 
221 
222 
223 #endif