EIC Software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
HelixHough.h
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file HelixHough.h
1 #ifndef __HELIX_HOUGH_H__
2 #define __HELIX_HOUGH_H__
3 
4 #include "fastvec.h"
5 #include "HelixRange.h"
6 #include "HelixResolution.h"
7 #include "SimpleHit3D.h"
8 #include "SimpleTrack3D.h"
9 #include "HelixKalmanState.h"
10 #include <xmmintrin.h>
11 #include <emmintrin.h>
12 #include <cmath>
13 #include <vector>
14 #include <set>
15 #include <map>
16 
17 #ifndef M_PI
18 #define M_PI 3.14159265358979323846
19 #endif
20 
21 
22 
23 
24 
25 class ParRange
26 {
27 public:
28  ParRange(){}
29  ParRange(unsigned int min_cAngle, unsigned int max_cAngle, unsigned int min_dca, unsigned int max_dca, unsigned int min_invR, unsigned int max_invR, unsigned int min_theta, unsigned int max_theta, unsigned int min_z, unsigned int max_z) : min_k(min_invR), max_k(max_invR), min_phi(min_cAngle), max_phi(max_cAngle), min_d(min_dca), max_d(max_dca), min_dzdl(min_theta), max_dzdl(max_theta), min_z0(min_z), max_z0(max_z) {}
30  ~ParRange() {}
31 
32  void mergeRange(unsigned int phi, unsigned int d, unsigned int k, unsigned int dzdl, unsigned int z0)
33  {
34  min_phi = (((phi < min_phi)-1)&min_phi) ^ (((phi >= min_phi)-1)&phi);
35  max_phi = (((phi > max_phi)-1)&max_phi) ^ (((phi <= max_phi)-1)&phi);
36  min_d = (((d < min_d)-1)&min_d) ^ (((d >= min_d)-1)&d);
37  max_d = (((d > max_d)-1)&max_d) ^ (((d <= max_d)-1)&d);
38  min_k = (((k < min_k)-1)&min_k) ^ (((k >= min_k)-1)&k);
39  max_k = (((k > max_k)-1)&max_k) ^ (((k <= max_k)-1)&k);
40  min_dzdl = (((dzdl < min_dzdl)-1)&min_dzdl) ^ (((dzdl >= min_dzdl)-1)&dzdl);
41  max_dzdl = (((dzdl > max_dzdl)-1)&max_dzdl) ^ (((dzdl <= max_dzdl)-1)&dzdl);
42  min_z0 = (((z0 < min_z0)-1)&min_z0) ^ (((z0 >= min_z0)-1)&z0);
43  max_z0 = (((z0 > max_z0)-1)&max_z0) ^ (((z0 <= max_z0)-1)&z0);
44  }
45 
46  void mergeRange(ParRange const& other)
47  {
48  min_phi = (((other.min_phi < min_phi)-1)&min_phi) ^ (((other.min_phi >= min_phi)-1)&other.min_phi);
49  max_phi = (((other.max_phi > max_phi)-1)&max_phi) ^ (((other.max_phi <= max_phi)-1)&other.max_phi);
50  min_d = (((other.min_d < min_d)-1)&min_d) ^ (((other.min_d >= min_d)-1)&other.min_d);
51  max_d = (((other.max_d > max_d)-1)&max_d) ^ (((other.max_d <= max_d)-1)&other.max_d);
52  min_k = (((other.min_k < min_k)-1)&min_k) ^ (((other.min_k >= min_k)-1)&other.min_k);
53  max_k = (((other.max_k > max_k)-1)&max_k) ^ (((other.max_k <= max_k)-1)&other.max_k);
54  min_dzdl = (((other.min_dzdl < min_dzdl)-1)&min_dzdl) ^ (((other.min_dzdl >= min_dzdl)-1)&other.min_dzdl);
55  max_dzdl = (((other.max_dzdl > max_dzdl)-1)&max_dzdl) ^ (((other.max_dzdl <= max_dzdl)-1)&other.max_dzdl);
56  min_z0 = (((other.min_z0 < min_z0)-1)&min_z0) ^ (((other.min_z0 >= min_z0)-1)&other.min_z0);
57  max_z0 = (((other.max_z0 > max_z0)-1)&max_z0) ^ (((other.max_z0 <= max_z0)-1)&other.max_z0);
58  }
59 
60  unsigned int min_k = UINT_MAX, max_k = UINT_MAX;
61  unsigned int min_phi = UINT_MAX, max_phi = UINT_MAX;
62  unsigned int min_d = UINT_MAX, max_d = UINT_MAX;
63  unsigned int min_dzdl = UINT_MAX, max_dzdl = UINT_MAX;
64  unsigned int min_z0 = UINT_MAX, max_z0 = UINT_MAX;
65 };
66 
67 
68 class ParameterCluster
69 {
70  public:
71  ParameterCluster(){}
72  ~ParameterCluster(){}
73 
74  ParRange range;
75  std::vector<unsigned int> hit_indexes;
76 };
77 
78 class BinEntryPair5D
79 {
80 public:
81  BinEntryPair5D(unsigned int b=0, unsigned int e=0, bool is=false) : bin(b), entry(e), is_seed(is) {}
82  ~BinEntryPair5D(){}
83 
84  bool operator<(const BinEntryPair5D& other) const
85  {
86  return ( bin < other.bin );
87  }
88 
89  static unsigned int linearBin(unsigned int nbins2, unsigned int nbins3, unsigned int nbins4, unsigned int nbins5, unsigned int bin1, unsigned int bin2, unsigned int bin3, unsigned int bin4, unsigned int bin5)
90  {
91  return bin5 + nbins5*(bin4 + nbins4*(bin3 + nbins3*( bin2 + nbins2*bin1 )));
92  }
93 
94  void bin5D(unsigned int nbins2, unsigned int nbins3, unsigned int nbins4, unsigned int nbins5, unsigned int& bin1, unsigned int& bin2, unsigned int& bin3, unsigned int& bin4, unsigned int& bin5) const
95  {
96  unsigned int temp1 = nbins2*nbins3*nbins4*nbins5;
97  bin1 = bin/temp1;
98  unsigned int temp2 = bin - bin1*temp1;
99  temp1 = nbins3*nbins4*nbins5;
100  bin2 = temp2/temp1;
101  temp2 -= bin2*temp1;
102  temp1 = nbins4*nbins5;
103  bin3 = temp2/temp1;
104  temp2 -= bin3*temp1;
105  temp1 = nbins5;
106  bin4 = temp2/temp1;
107  temp2 -= bin4*temp1;
108  bin5 = temp2;
109  }
110 
111  unsigned int bin;
112  unsigned int entry;
113  bool is_seed;
114 };
115 
116 
117 
118 class HelixHough
119 {
120  public:
121  HelixHough(unsigned int n_phi, unsigned int n_d, unsigned int n_k, unsigned int n_dzdl, unsigned int n_z0, HelixResolution& min_resolution, HelixResolution& max_resolution, HelixRange& range);
122  HelixHough(std::vector<std::vector<unsigned int> >& zoom_profile, unsigned int minzoom, HelixRange& range);
123  virtual ~HelixHough();
124 
125  void initHelixHough(unsigned int n_phi, unsigned int n_d, unsigned int n_k, unsigned int n_dzdl, unsigned int n_z0, HelixResolution& min_resolution, HelixResolution& max_resolution, HelixRange& range);
126 
127  void findHelices(std::vector<SimpleHit3D>& hits, unsigned int min_hits, unsigned int max_hits, std::vector<SimpleTrack3D>& tracks, unsigned int maxtracks=0);
128 
129  void findHelices(unsigned int min_hits, unsigned int max_hits, std::vector<SimpleTrack3D>& tracks, unsigned int maxtracks, unsigned int zoomlevel);
130 
131  void findSeededHelices(std::vector<SimpleTrack3D>& seeds, std::vector<SimpleHit3D>& hits, unsigned int min_hits, unsigned int max_hits, std::vector<SimpleTrack3D>& tracks, unsigned int maxtracks=0);
132  void findSeededHelices_run(std::vector<SimpleTrack3D>& seeds, std::vector<SimpleHit3D>& hits, unsigned int min_hits, unsigned int max_hits, std::vector<SimpleTrack3D>& tracks, unsigned int maxtracks=0);
133  void findSeededHelices(unsigned int min_hits, unsigned int max_hits, std::vector<SimpleTrack3D>& tracks, unsigned int maxtracks, unsigned int zoomlevel);
134 
135  void vote(unsigned int zoomlevel);
136  static void phiRange_sse(float* hit_x, float* hit_y, float* min_d, float* max_d, float* min_k, float* max_k, float* min_phi_1, float* max_phi_1, float* min_phi_2, float* max_phi_2);
137  static void phiRange_sse(float* hit_x, float* hit_y, float* min_d, float* max_d, float* min_k, float* max_k, float* min_phi, float* max_phi, float hel, __m128& phi_3_out, __m128& phi_4_out);
138  static void phiRange_sse(float* hit_x, float* hit_y, float* min_d, float* max_d, float* max_k, float* min_phi, float* max_phi, float hel, __m128& phi_3, __m128& phi_4, __m128& phi_3_out, __m128& phi_4_out);
139  static void dzdlRange_sse(float* x_a, float* y_a, float* z_a, float cosphi1, float sinphi1, float cosphi2, float sinphi2, float min_k, float max_k, float min_d, float max_d, float* min_z0, float* max_z0, float* min_dzdl_a, float* max_dzdl_a);
140  static void phiRange_sse(float* hit_x, float* hit_y, float* min_d, float* max_d, float* min_k, float* max_k, float* min_phi, float* max_phi, float* min_phi_2, float* max_phi_2, float hel, __m128& phi_3_out, __m128& phi_4_out, float* hit_x_2, float* hit_y_2, __m128& phi_3_out_2, __m128& phi_4_out_2);
141  static void phiRange_sse(float* hit_x, float* hit_y, float* min_d, float* max_d, float* min_k, float* max_k, float* min_phi, float* max_phi, float* min_phi_2, float* max_phi_2, float hel, __m128& phi_3, __m128& phi_4, __m128& phi_3_out, __m128& phi_4_out, float* hit_x_2, float* hit_y_2, __m128& phi_3_2, __m128& phi_4_2, __m128& phi_3_out_2, __m128& phi_4_out_2);
142 
143  void setPrintTimings(bool pt){print_timings=pt;}
144 
145  virtual void finalize(std::vector<SimpleTrack3D>&, std::vector<SimpleTrack3D>&){}
146  virtual void findTracks(std::vector<SimpleHit3D>&, std::vector<SimpleTrack3D>&, const HelixRange&) = 0;
147  virtual void initEvent(std::vector<SimpleHit3D>&, unsigned int){}
148  virtual void findSeededTracks(std::vector<SimpleTrack3D>&, std::vector<SimpleHit3D>&, std::vector<SimpleTrack3D>&, const HelixRange&){}
149 
150  // return true if we want to go straight into the user-written findTracks function instead of possibly continuing the Hough Transform
151  virtual bool breakRecursion(const std::vector<SimpleHit3D>&, const HelixRange&){return false;}
152 
153  // additional phi error to add to a hit in the voting stage, given an bin with the k and dzdl parameters passed. This is useful
154  // to take into account multiple scattering
155  virtual float phiError(SimpleHit3D&, float, float, float, float, float, float, float, float, bool = false){return 0.;}
156  // same as above, but for additional error on dzdl
157  virtual float dzdlError(SimpleHit3D&, float, float, float, float, float, float, float, float, bool = false){return 0.;}
158 
159  virtual void initSeeding(std::vector<SimpleTrack3D>&){}
160 
161  virtual void setSeparateByHelicity(bool sbh){separate_by_helicity=sbh;}
162  virtual void setOnlyOneHelicity(bool ooh){only_one_helicity=ooh;}
163  void setHelicity(bool hel){helicity=hel;}
164 
165  virtual void requireLayers(unsigned int nl)
166  {
167  check_layers = true;
168  req_layers = nl;
169  }
170 
171  virtual void setBinScale(float b_scl){bin_scale = b_scl;}
172  virtual void setZBinScale(float b_scl){z_bin_scale = b_scl;}
173  virtual void setRemoveHits(bool rh){remove_hits=rh;}
174 
175  virtual void setRangeFromSeed(HelixRange& range, SimpleTrack3D& seed);
176 
177  virtual void setTopRange(HelixRange& tr){top_range = tr;}
178  void splitIntoBins(unsigned int min_hits, unsigned int max_hits, std::vector<HelixRange>& ranges, std::vector<std::vector<SimpleHit3D> >& split_hits, unsigned int zoomlevel);
179 
180  virtual void clear(){}
181 
182  virtual void setStartZoom(unsigned int sz){start_zoom=sz;}
183 
184  virtual void setMaxHitsPairs(unsigned int mhp){max_hits_pairs=mhp;}
185 
186  void setClusterStartBin(unsigned int csb){cluster_start_bin=csb;}
187 
188  void setSeedStates(std::vector<HelixKalmanState>& states){seed_states = states;}
189 
190  std::vector<HelixKalmanState>& getKalmanStates(){return track_states;}
191 
192  void setLayersAtATime(unsigned int l){layers_at_a_time = l;}
193 
194  void setSmoothBack(bool sb){smooth_back=sb;}
195 
196  void setCullInputHits( bool cih ){ cull_input_hits = cih; }
197  void setIterateClustering( bool icl ){ iterate_clustering = icl; }
198 
199  protected:
200  bool remove_hits;
201  std::vector<unsigned int>* hit_used;
202 
203  void setRange(const BinEntryPair5D& bp, HelixRange& range1, HelixRange& range2, unsigned int n_phi, unsigned int n_d, unsigned int n_k, unsigned int n_dzdl, unsigned int n_z0);
204 
205  //number of bins for each zoom level
206  std::vector<unsigned int> n_phi_bins;
207  std::vector<unsigned int> n_d_bins;
208  std::vector<unsigned int> n_k_bins;
209  std::vector<unsigned int> n_dzdl_bins;
210  std::vector<unsigned int> n_z0_bins;
211  //top level hits
212  std::vector<SimpleHit3D>* base_hits;
213  // vector of hits used for isolating hits in a single zoom level
214  std::vector<std::vector<SimpleHit3D>* > hits_vec;
215  // for each zoomlevel, a vector of pair indexes, with the index mapping to the entry position in hits_vec
216  std::vector<std::vector<std::pair<unsigned int,unsigned int> >* > pairs_vec;
217  // stores external index values
218  std::vector<unsigned int> index_mapping;
220  std::vector<std::vector<BinEntryPair5D>* > bins_vec;
221 
222  std::vector<std::vector<SimpleTrack3D>* > seeds_vec;
223 
224  HelixRange top_range;
225  HelixRange current_range;
226  std::vector<HelixRange> zoomranges;
227  std::vector<std::vector<ParameterCluster>*> clusters_vec;
228  std::vector<unsigned int> num_clusters;
229  // vector on which to perform unique sorting
230  std::vector<unsigned int> C_clus;
231  // vector to contain merged hits from a cluster and candidate bin
232  std::vector<unsigned int> temp_merged_clus;
233 
234  unsigned int max_zoom;//maximum number of times to zoom in on the parameter array
235  unsigned int min_zoom;//minimum number of times to zoom in on the parameter array
236 
237  bool using_vertex;
238 
239  unsigned int max_tracks;
240 
241  double vote_time;
242  double xy_vote_time;
243  double z_vote_time;
244  double cluster_time;
245 
246  bool print_timings;
247 
248  bool separate_by_helicity,helicity,only_one_helicity;
249 
250  bool check_layers;
251  unsigned int req_layers;
252 
253  float bin_scale; // proportion of the bin size in each direction which is actually used for voting
254  float z_bin_scale; // proportion of the bin size in the z-direction in each direction which is actually used for voting
255 
256  unsigned int start_zoom; // top-level zoomlevel : defaults to zero
257 
258 
259  static void allButKappaRange_sse(float* x1_a,float* x2_a,float* y1_a,float* y2_a,float* z1_a,float* z2_a, float* min_k_a,float* max_k_a, float* min_phi_1_a,float* max_phi_1_a,float* min_phi_2_a,float* max_phi_2_a, float* min_d_1_a,float* max_d_1_a,float* min_d_2_a,float* max_d_2_a, float* min_dzdl_a,float* max_dzdl_a, float* min_z0_1_a,float* max_z0_1_a,float* min_z0_2_a,float* max_z0_2_a);
260 
261 
262  void fillBins(unsigned int total_bins, unsigned int hit_counter, float* min_phi_a, float* max_phi_a, std::vector<SimpleHit3D>& four_hits, fastvec2d& z_bins, unsigned int n_d, unsigned int n_k, unsigned int n_dzdl, unsigned int n_z0, unsigned int d_bin, unsigned int k_bin, unsigned int n_phi, unsigned int zoomlevel, float low_phi, float high_phi, float inv_phi_range, fastvec& vote_array);
263 
264  void makeClusters(unsigned int zoomlevel, unsigned int MAX, unsigned int n_phi, unsigned int n_d, unsigned int n_k, unsigned int n_dzdl, unsigned int n_z0, unsigned int min_hits, std::vector<ParameterCluster>& clusters, bool& use_clusters, bool& is_super_bin);
265 
266  bool attemptClusterMerge(unsigned int zoomlevel, unsigned int MAX, unsigned int ca, unsigned int d, unsigned int r, unsigned int th, unsigned int zz0, unsigned int bin, unsigned int newbin, std::vector<unsigned char>& good_bins, unsigned int volume, float cluster_size_cut, float overlap_cut, std::vector<ParameterCluster>& clusters, unsigned int* bins_start, unsigned int* bins_end, std::vector<unsigned int>& map_clus, std::vector<unsigned char>& too_big, std::vector<unsigned int>& temp_merged, std::vector<unsigned int>& C);
267 
268  void vote_z(unsigned int zoomlevel, unsigned int n_phi, unsigned int n_d, unsigned int
269  n_k, unsigned int n_dzdl, unsigned int n_z0, fastvec2d& z_bins);
270 
271 
272  void vote_pairs(unsigned int zoomlevel);
273  void fillBins(unsigned int total_bins, unsigned int pair_counter, unsigned int* pair_index, float* min_phi, float* max_phi, float* min_d, float* max_d, float* min_dzdl, float* max_dzdl, float* min_z0, float* max_z0, std::vector<std::vector<SimpleHit3D> > & four_pairs, unsigned int n_d, unsigned int n_k, unsigned int n_dzdl, unsigned int n_z0, unsigned int k_bin, unsigned int n_phi, unsigned int zoomlevel, float low_phi, float high_phi, float low_d, float high_d, float low_z0, float high_z0, float low_dzdl, float high_dzdl, float inv_phi_range, float inv_d_range, float inv_z0_range, float inv_dzdl_range, fastvec& vote_array);
274  void findHelicesByPairsBegin(unsigned int min_hits, unsigned int max_hits, std::vector<SimpleTrack3D>& tracks, unsigned int maxtracks, unsigned int zoomlevel);
275  void findHelicesByPairs(unsigned int min_hits, unsigned int max_hits, std::vector<SimpleTrack3D>& tracks, unsigned int maxtracks, unsigned int zoomlevel);
276 
277  unsigned int max_hits_pairs;
278 
279  std::vector<std::pair<unsigned int, unsigned int> > temp_pairs;
280  std::set<unsigned int> new_hits;
281  std::map<unsigned int, unsigned int> old_to_new;
282 
283  unsigned int cluster_start_bin;
284 
285  unsigned int bins_start[1<<12];
286  unsigned int bins_end[1<<12];
287 
288  unsigned int layers_at_a_time;
289 
290  std::vector<HelixKalmanState> track_states;
291  std::vector<HelixKalmanState> seed_states;
292 
293  unsigned int n_layers;
294  int layer_start;
295  int layer_end;
296  bool smooth_back;
297  bool cull_input_hits;
298  bool iterate_clustering;
299 };
300 
301 #endif
302