EIC Software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
HelixHough_findHelices.cpp
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file HelixHough_findHelices.cpp
1 #include "HelixHough.h"
2 #include "HelixRange.h"
3 #include "SimpleHit3D.h"
4 #include "SimpleTrack3D.h"
5 
6 #include <algorithm>
7 #include <cstddef>
8 #include <cmath>
9 #include <iostream>
10 #include <memory>
11 #include <sys/time.h>
12 #include <utility>
13 #include <vector>
14 
15 using namespace std;
16 
17 // static inline void in_place_counting_sort_unique(vector<unsigned int>& A,
18 // vector<unsigned int>& C,
19 // unsigned int MAX) {
20 // unsigned int SIZE = A.size();
21 // for (unsigned int i = 0; i < SIZE; ++i) {
22 // ++C[A[i]];
23 // }
24 // unsigned int current = 0;
25 // for (unsigned int i = 0; i < MAX; ++i) {
26 // A[current] = ((((C[i] != 0) - 1) & (A[current])) ^ (((C[i] == 0) - 1) & i));
27 // current += (C[i] != 0);
28 // }
29 // A.resize(current);
30 // }
31 
32 static inline void in_place_counting_unique(vector<unsigned int>& A,
33  vector<unsigned int>& C) {
34  unsigned int SIZE = A.size();
35  unsigned int current = 0;
36  for (unsigned int i = 0; i < SIZE; ++i) {
37  C[A[i]] += 1;
38  if (C[A[i]] == 1) {
39  A[current] = A[i];
40  current += 1;
41  }
42  }
43  A.resize(current);
44 }
45 
46 void HelixHough::setRange(const BinEntryPair5D& bp, HelixRange& range1,
47  HelixRange& range2, unsigned int n_phi,
48  unsigned int n_d, unsigned int n_k,
49  unsigned int n_dzdl, unsigned int n_z0) {
50  float dzdl_size = (range1.max_dzdl - range1.min_dzdl) / ((float)n_dzdl);
51  float z0_size = (range1.max_z0 - range1.min_z0) / ((float)n_z0);
52  float phi_size = (range1.max_phi - range1.min_phi) / ((float)n_phi);
53  float d_size = (range1.max_d - range1.min_d) / ((float)n_d);
54  float k_size = (range1.max_k - range1.min_k) / ((float)n_k);
55 
56  unsigned int z0_bin = 0;
57  unsigned int dzdl_bin = 0;
58  unsigned int k_bin = 0;
59  unsigned int d_bin = 0;
60  unsigned int phi_bin = 0;
61 
62  bp.bin5D(n_d, n_k, n_dzdl, n_z0, phi_bin, d_bin, k_bin, dzdl_bin, z0_bin);
63  range2.min_phi = range1.min_phi + phi_size * ((float)(phi_bin));
64  range2.max_phi = range2.min_phi + phi_size;
65  range2.min_d = range1.min_d + d_size * ((float)(d_bin));
66  range2.max_d = range2.min_d + d_size;
67  range2.min_k = range1.min_k + k_size * ((float)(k_bin));
68  range2.max_k = range2.min_k + k_size;
69  range2.min_dzdl = range1.min_dzdl + dzdl_size * ((float)(dzdl_bin));
70  range2.max_dzdl = range2.min_dzdl + dzdl_size;
71  range2.min_z0 = range1.min_z0 + z0_size * ((float)(z0_bin));
72  range2.max_z0 = range2.min_z0 + z0_size;
73 }
74 
75 static inline void setClusterRange(HelixRange& range1, HelixRange& range2,
76  ParRange& prange, unsigned int n_phi,
77  unsigned int n_d, unsigned int n_k,
78  unsigned int n_dzdl, unsigned int n_z0) {
79  float dzdl_size = (range1.max_dzdl - range1.min_dzdl) / ((float)n_dzdl);
80  float z0_size = (range1.max_z0 - range1.min_z0) / ((float)n_z0);
81  float phi_size = (range1.max_phi - range1.min_phi) / ((float)n_phi);
82  float d_size = (range1.max_d - range1.min_d) / ((float)n_d);
83  float k_size = (range1.max_k - range1.min_k) / ((float)n_k);
84 
85  range2.min_phi = range1.min_phi + phi_size * ((float)(prange.min_phi));
86  range2.max_phi = range1.min_phi + phi_size * ((float)(prange.max_phi + 1));
87  range2.min_d = range1.min_d + d_size * ((float)(prange.min_d));
88  range2.max_d = range1.min_d + d_size * ((float)(prange.max_d + 1));
89  range2.min_k = range1.min_k + k_size * ((float)(prange.min_k));
90  range2.max_k = range1.min_k + k_size * ((float)(prange.max_k + 1));
91  range2.min_dzdl = range1.min_dzdl + dzdl_size * ((float)(prange.min_dzdl));
92  range2.max_dzdl =
93  range1.min_dzdl + dzdl_size * ((float)(prange.max_dzdl + 1));
94  range2.min_z0 = range1.min_z0 + z0_size * ((float)(prange.min_z0));
95  range2.max_z0 = range1.min_z0 + z0_size * ((float)(prange.max_z0 + 1));
96 }
97 
98 void HelixHough::findHelices(vector<SimpleHit3D>& hits_init, unsigned int min_hits,
99  unsigned int max_hits,
100  vector<SimpleTrack3D>& tracks,
101  unsigned int maxtracks) {
102  vote_time = 0.;
103  xy_vote_time = 0.;
104  z_vote_time = 0.;
105  cluster_time = 0.;
106 
107  vector<SimpleHit3D> hits;
108 
109  if (layer_end > 0 && layer_start > 0){
110  for (unsigned int in = 0; in < hits_init.size(); in++) {
111  if(hits_init[in].get_layer() > layer_start-1 && hits_init[in].get_layer() < layer_end+1 )
112  hits.push_back(hits_init[in]);
113  }
114  } else if ( layer_end < 0 && layer_start > 0) {
115  for (unsigned int in = 0; in < hits_init.size(); in++) {
116  if(hits_init[in].get_layer() > layer_start-1)
117  hits.push_back(hits_init[in]);
118  }
119  } else if ( layer_end > 0 && layer_start < 0) {
120  for (unsigned int in = 0; in < hits_init.size(); in++) {
121  if( hits_init[in].get_layer() < layer_end+1 )
122  hits.push_back(hits_init[in]);
123  }
124  } else {
125  for (unsigned int in = 0; in < hits_init.size(); in++) {
126  hits.push_back(hits_init[in]);
127  }
128  }
129 
130  index_mapping.clear();
131  index_mapping.resize(hits.size(), 0);
132  hit_used->clear();
133  for (unsigned int i = 0; i < hits.size(); i++) {
134  index_mapping[i] = hits[i].get_id();
135  hits[i].set_id(i);
136  }
137  (*hit_used).assign(hits.size(), false);
138 
139  initEvent(hits, min_hits);
140 
141  max_tracks = maxtracks;
142  base_hits = &hits;
143  (*(hits_vec[start_zoom])) = hits;
144  current_range = top_range;
145  zoomranges.clear();
146  for (unsigned int z = 0; z <= max_zoom; z++) {
147  zoomranges.push_back(top_range);
148  }
149  vector<SimpleTrack3D> temp_tracks;
150  if ((separate_by_helicity == true) && (only_one_helicity == false)) {
151  helicity = true;
152  findHelices(min_hits, max_hits, temp_tracks, maxtracks, start_zoom);
153  helicity = false;
154  findHelices(min_hits, max_hits, temp_tracks, maxtracks, start_zoom);
155  } else {
156  findHelices(min_hits, max_hits, temp_tracks, maxtracks, start_zoom);
157  }
158 
159  vector<SimpleHit3D> tr_hits;
160  if (cull_input_hits == true) {
161  for (unsigned int i = 0; i < hits.size(); i++) {
162  if ((*hit_used)[hits[i].get_id()] == false) {
163  tr_hits.push_back(hits[i]);
164  tr_hits.back().set_id(index_mapping[i]);
165  }
166  }
167  }
168 
169  for (unsigned int i = 0; i < hits.size(); i++) {
170  hits[i].set_id(index_mapping[i]);
171  }
172  for (unsigned int t = 0; t < temp_tracks.size(); t++) {
173  for (unsigned int h = 0; h < temp_tracks[t].hits.size(); h++) {
174  if (temp_tracks[t].hits[h].get_id() != (unsigned)-1) {
175  temp_tracks[t].hits[h].set_id(index_mapping[temp_tracks[t].hits[h].get_id()]);
176  }
177  }
178  }
179 
180  finalize(temp_tracks, tracks);
181 
182  if (cull_input_hits == true) {
183  hits = tr_hits;
184  }
185 
186  if (print_timings == true) {
187  cout << "vote time = " << vote_time << endl;
188  cout << "xy vote time = " << xy_vote_time << endl;
189  cout << "z vote time = " << z_vote_time << endl;
190  cout << "cluster time = " << cluster_time << endl;
191  }
192 }
193 
194 static bool mergeClusters(ParameterCluster& clus1,
195  ParameterCluster const& clus2, unsigned int MAX,
196  float overlap_cut) {
197  vector<unsigned int> old_indexes = clus1.hit_indexes;
198  for (unsigned int i = 0; i < clus2.hit_indexes.size(); ++i) {
199  clus1.hit_indexes.push_back(clus2.hit_indexes[i]);
200  }
201  vector<unsigned int> C;
202  C.assign(MAX + 1, 0);
203  in_place_counting_unique(clus1.hit_indexes, C);
204 
205  float size_diff_1 = ((float)(clus1.hit_indexes.size() - old_indexes.size()) /
206  ((float)(old_indexes.size())));
207  float size_diff_2 =
208  ((float)(clus1.hit_indexes.size() - clus2.hit_indexes.size()) /
209  ((float)(clus2.hit_indexes.size())));
210 
211  if ((size_diff_1 < overlap_cut) || (size_diff_2 < overlap_cut)) {
212  clus1.range.mergeRange(clus2.range);
213  return true;
214  } else {
215  clus1.hit_indexes = old_indexes;
216  return false;
217  }
218 }
219 
220 bool HelixHough::attemptClusterMerge(
221  unsigned int zoomlevel, unsigned int MAX, unsigned int ca, unsigned int d,
222  unsigned int r, unsigned int th, unsigned int zz0, unsigned int bin,
223  unsigned int newbin, vector<unsigned char>& good_bins, unsigned int volume,
224  float cluster_size_cut, float overlap_cut,
225  vector<ParameterCluster>& clusters, unsigned int* bins_start,
226  unsigned int* bins_end, vector<unsigned int>& map_clus,
227  vector<unsigned char>& too_big, vector<unsigned int>& /*temp_merged*/,
228  vector<unsigned int>& C) {
229  if ((good_bins[newbin] == 1) && (map_clus[newbin] != 4294967295)) {
230  if (too_big[map_clus[newbin]] == 0) {
231  unsigned int tempsize = clusters[map_clus[newbin]].hit_indexes.size();
232  for (unsigned int ind = bins_start[bin]; ind <= bins_end[bin]; ++ind) {
233  clusters[map_clus[newbin]].hit_indexes.push_back(
234  (*(bins_vec[zoomlevel]))[ind].entry);
235  }
236  C.clear();
237  C.assign(MAX + 1, 0);
238  in_place_counting_unique(clusters[map_clus[newbin]].hit_indexes, C);
239  unsigned int size_diff =
240  clusters[map_clus[newbin]].hit_indexes.size() - tempsize;
241  unsigned int size_diff_2 = clusters[map_clus[newbin]].hit_indexes.size() -
242  (1 + bins_end[bin] - bins_start[bin]);
243  if ((((float)size_diff) / ((float)(1 + bins_end[bin] - bins_start[bin])) <
244  overlap_cut) ||
245  (((float)size_diff_2) / ((float)(tempsize)) < overlap_cut)) {
246  clusters[map_clus[newbin]].range.mergeRange(ca, d, r, th, zz0);
247  ParameterCluster* cluster = &(clusters[map_clus[newbin]]);
248  unsigned int cluster_volume =
249  ((cluster->range.max_phi - cluster->range.min_phi) + 1) *
250  ((cluster->range.max_d - cluster->range.min_d) + 1) *
251  ((cluster->range.max_k - cluster->range.min_k) + 1) *
252  ((cluster->range.max_dzdl - cluster->range.min_dzdl) + 1) *
253  ((cluster->range.max_z0 - cluster->range.min_z0) + 1);
254  if (((float)cluster_volume) / ((float)volume) > cluster_size_cut) {
255  too_big[map_clus[newbin]] = 1;
256  }
257  map_clus[bin] = map_clus[newbin];
258  return true;
259  } else {
260  clusters[map_clus[newbin]].hit_indexes.resize(tempsize);
261  }
262  }
263  }
264  return false;
265 }
266 
267 void HelixHough::makeClusters(unsigned int zoomlevel, unsigned int MAX,
268  unsigned int n_phi, unsigned int n_d,
269  unsigned int n_k, unsigned int n_dzdl,
270  unsigned int n_z0, unsigned int min_hits,
271  vector<ParameterCluster>& clusters,
272  bool& /*use_clusters*/, bool& is_super_bin) {
273  unsigned int volume = n_phi * n_d * n_k * n_dzdl * n_z0;
274  float cluster_size_cut = 1.0;
275 // float bin_size_cut = 0.75;
276  float overlap_cut = 0.1;
277  is_super_bin = false;
278 
279  vector<unsigned int> map_clus(volume, 4294967295);
280  vector<unsigned char> good_bins(volume, 0);
281  vector<unsigned char> too_big(volume, 0);
282 
283  for (unsigned int ca = 0; ca < n_phi; ++ca) {
284  for (unsigned int d = 0; d < n_d; ++d) {
285  for (unsigned int r = 0; r < n_k; ++r) {
286  for (unsigned int th = 0; th < n_dzdl; ++th) {
287  for (unsigned int zz0 = 0; zz0 < n_z0; ++zz0) {
288  unsigned int bin = BinEntryPair5D::linearBin(n_d, n_k, n_dzdl, n_z0,
289  ca, d, r, th, zz0);
290 
291  if (bins_end[bin] == 4294967295) {
292  continue;
293  }
294  if ((1 + bins_end[bin] - bins_start[bin]) >= min_hits) {
295  if (check_layers == true) {
296  unsigned int layer_mask[4] = {0, 0, 0, 0};
297  for (unsigned int i = bins_start[bin]; i <= bins_end[bin];
298  ++i) {
299  if ((*(hits_vec[zoomlevel]))[(*(bins_vec[zoomlevel]))[i]
300  .entry]
301  .get_layer() < 32) {
302  layer_mask[0] =
303  layer_mask[0] |
304  (1 << (*(
305  hits_vec[zoomlevel]))[(*(bins_vec[zoomlevel]))[i]
306  .entry]
307  .get_layer());
308  } else if ((*(hits_vec[zoomlevel]))
309  [(*(bins_vec[zoomlevel]))[i].entry]
310  .get_layer() < 64) {
311  layer_mask[1] =
312  layer_mask[1] |
313  (1 << ((*(hits_vec[zoomlevel]))
314  [(*(bins_vec[zoomlevel]))[i].entry]
315  .get_layer() -
316  32));
317  } else if ((*(hits_vec[zoomlevel]))
318  [(*(bins_vec[zoomlevel]))[i].entry]
319  .get_layer() < 96) {
320  layer_mask[2] =
321  layer_mask[2] |
322  (1 << ((*(hits_vec[zoomlevel]))
323  [(*(bins_vec[zoomlevel]))[i].entry]
324  .get_layer() -
325  64));
326  } else if ((*(hits_vec[zoomlevel]))
327  [(*(bins_vec[zoomlevel]))[i].entry]
328  .get_layer() < 128) {
329  layer_mask[3] =
330  layer_mask[3] |
331  (1 << ((*(hits_vec[zoomlevel]))
332  [(*(bins_vec[zoomlevel]))[i].entry]
333  .get_layer() -
334  96));
335  }
336  }
337  unsigned int nlayers = __builtin_popcount(layer_mask[0]) +
338  __builtin_popcount(layer_mask[1]) +
339  __builtin_popcount(layer_mask[2]) +
340  __builtin_popcount(layer_mask[3]);
341  if (nlayers >= req_layers) {
342  good_bins[bin] = 1;
343  }
344  } else {
345  good_bins[bin] = 1;
346  }
347  } else {
348  continue;
349  }
350 
351  if (good_bins[bin] == 0) {
352  continue;
353  }
354 
355  if (ca > 0) {
356  unsigned int newbin = bin - n_d * n_k * n_dzdl * n_z0;
357  if (attemptClusterMerge(
358  zoomlevel, MAX, ca, d, r, th, zz0, bin, newbin, good_bins,
359  volume, cluster_size_cut, overlap_cut, clusters,
360  (unsigned int*)(bins_start), (unsigned int*)(bins_end),
361  map_clus, too_big, temp_merged_clus, C_clus) == true) {
362  continue;
363  }
364  }
365 
366  if (d > 0) {
367  unsigned int newbin = bin - n_k * n_dzdl * n_z0;
368  if (attemptClusterMerge(
369  zoomlevel, MAX, ca, d, r, th, zz0, bin, newbin, good_bins,
370  volume, cluster_size_cut, overlap_cut, clusters,
371  (unsigned int*)(bins_start), (unsigned int*)(bins_end),
372  map_clus, too_big, temp_merged_clus, C_clus) == true) {
373  continue;
374  }
375  }
376 
377  if (r > 0) {
378  unsigned int newbin = bin - n_dzdl * n_z0;
379  if (attemptClusterMerge(
380  zoomlevel, MAX, ca, d, r, th, zz0, bin, newbin, good_bins,
381  volume, cluster_size_cut, overlap_cut, clusters,
382  (unsigned int*)(bins_start), (unsigned int*)(bins_end),
383  map_clus, too_big, temp_merged_clus, C_clus) == true) {
384  continue;
385  }
386  }
387 
388  if (th > 0) {
389  unsigned int newbin = bin - n_z0;
390  if (attemptClusterMerge(
391  zoomlevel, MAX, ca, d, r, th, zz0, bin, newbin, good_bins,
392  volume, cluster_size_cut, overlap_cut, clusters,
393  (unsigned int*)(bins_start), (unsigned int*)(bins_end),
394  map_clus, too_big, temp_merged_clus, C_clus) == true) {
395  continue;
396  }
397  }
398 
399  if (zz0 > 0) {
400  unsigned int newbin = bin - 1;
401  if (attemptClusterMerge(
402  zoomlevel, MAX, ca, d, r, th, zz0, bin, newbin, good_bins,
403  volume, cluster_size_cut, overlap_cut, clusters,
404  (unsigned int*)(bins_start), (unsigned int*)(bins_end),
405  map_clus, too_big, temp_merged_clus, C_clus) == true) {
406  continue;
407  }
408  }
409  if (num_clusters[zoomlevel] >= clusters.size()) {
410  clusters.push_back(ParameterCluster());
411  }
412  clusters[num_clusters[zoomlevel]].range =
413  ParRange(ca, ca, d, d, r, r, th, th, zz0, zz0);
414  map_clus[bin] = num_clusters[zoomlevel];
415  clusters[num_clusters[zoomlevel]].hit_indexes.clear();
416  for (unsigned int ind = bins_start[bin]; ind <= bins_end[bin];
417  ++ind) {
418  clusters[num_clusters[zoomlevel]].hit_indexes.push_back(
419  (*(bins_vec[zoomlevel]))[ind].entry);
420  }
421 
422  num_clusters[zoomlevel] += 1;
423  }
424  }
425  }
426  }
427  }
428  if (iterate_clustering == false) {
429  return;
430  }
431  if (num_clusters[zoomlevel] == 0) {
432  return;
433  }
434  vector<ParameterCluster> in_clusters;
435  for (unsigned int i = 0; i < num_clusters[zoomlevel]; ++i) {
436  in_clusters.push_back(clusters[i]);
437  }
438  vector<ParameterCluster> out_clusters;
439  out_clusters.push_back(in_clusters[0]);
440  for (unsigned int i = 1; i < in_clusters.size(); ++i) {
441  bool merged = false;
442  for (unsigned int j = 0; j < out_clusters.size(); ++j) {
443  merged = mergeClusters(out_clusters[j], in_clusters[i], MAX, overlap_cut);
444  if (merged == true) {
445  merged = true;
446  break;
447  }
448  }
449  if (merged == false) {
450  out_clusters.push_back(in_clusters[i]);
451  }
452  }
453  clusters = out_clusters;
454  num_clusters[zoomlevel] = out_clusters.size();
455 }
456 
457 void HelixHough::findHelices(unsigned int min_hits, unsigned int max_hits,
458  vector<SimpleTrack3D>& tracks,
459  unsigned int maxtracks, unsigned int zoomlevel) {
460  unsigned int tracks_at_start = tracks.size();
461 
462  if ((maxtracks != 0) && (tracks.size() >= max_tracks)) {
463  return;
464  }
465 
466  timeval t1, t2;
467  double time1 = 0.;
468  double time2 = 0.;
469  if (print_timings == true) {
470  gettimeofday(&t1, NULL);
471  }
472  vote(zoomlevel);
473  if (print_timings == true) {
474  gettimeofday(&t2, NULL);
475  time1 = ((double)(t1.tv_sec) + (double)(t1.tv_usec) / 1000000.);
476  time2 = ((double)(t2.tv_sec) + (double)(t2.tv_usec) / 1000000.);
477  vote_time += (time2 - time1);
478  }
479 
480  unsigned int n_entries = bins_vec[zoomlevel]->size();
481  if (n_entries == 0) {
482  return;
483  }
484 
485  unsigned int n_phi = n_phi_bins[zoomlevel];
486  unsigned int n_d = n_d_bins[zoomlevel];
487  unsigned int n_k = n_k_bins[zoomlevel];
488  unsigned int n_dzdl = n_dzdl_bins[zoomlevel];
489  unsigned int n_z0 = n_z0_bins[zoomlevel];
490 
491  num_clusters[zoomlevel] = 0;
492  bool use_clusters = true;
493  bool is_super_bin = false;
494  if (zoomlevel >= cluster_start_bin) {
495  if (print_timings == true) {
496  gettimeofday(&t1, NULL);
497  }
498  makeClusters(zoomlevel, hits_vec[zoomlevel]->size(), n_phi, n_d, n_k,
499  n_dzdl, n_z0, min_hits, *(clusters_vec[zoomlevel]),
500  use_clusters, is_super_bin);
501 
502  if (print_timings) {
503  gettimeofday(&t2, NULL);
504  time1 = ((double)(t1.tv_sec) + (double)(t1.tv_usec) / 1000000.);
505  time2 = ((double)(t2.tv_sec) + (double)(t2.tv_usec) / 1000000.);
506  cluster_time += (time2 - time1);
507  }
508  } else {
509  use_clusters = false;
510  }
511 
512  if (use_clusters == false) {
513  unsigned int count = 0;
514  hits_vec[zoomlevel + 1]->clear();
515  setRange((*(bins_vec[zoomlevel]))[count], zoomranges[zoomlevel],
516  zoomranges[zoomlevel + 1], n_phi, n_d, n_k, n_dzdl, n_z0);
517  // scan over the bins in 5-D hough space
518  while (count < n_entries) {
519  if (remove_hits == true) {
520  if ((*hit_used)[(*(
521  hits_vec[zoomlevel]))[(*(bins_vec[zoomlevel]))[count].entry]
522  .get_id()] == false) {
523  hits_vec[zoomlevel + 1]->push_back(
524  (*(hits_vec[zoomlevel]))[(*(bins_vec[zoomlevel]))[count].entry]);
525  }
526  } else {
527  hits_vec[zoomlevel + 1]->push_back(
528  (*(hits_vec[zoomlevel]))[(*(bins_vec[zoomlevel]))[count].entry]);
529  }
530 
531  count += 1;
532  // we have collected all hits from this bin. now zoom again or find tracks
533  // with user routine
534  if ((count == n_entries) || ((*(bins_vec[zoomlevel]))[count].bin !=
535  (*(bins_vec[zoomlevel]))[count - 1].bin)) {
536  if (hits_vec[zoomlevel + 1]->size() >= min_hits) {
537  if (breakRecursion(*(hits_vec[zoomlevel + 1]),
538  zoomranges[zoomlevel + 1]) == true) {
539  } else if (((zoomlevel + 1) == max_zoom)) {
540  findTracks(*(hits_vec[zoomlevel + 1]), tracks,
541  zoomranges[zoomlevel + 1]);
542  } else if (((zoomlevel + 1) >= min_zoom) &&
543  ((hits_vec[zoomlevel + 1]->size() <= max_hits))) {
544  findTracks(*(hits_vec[zoomlevel + 1]), tracks,
545  zoomranges[zoomlevel + 1]);
546  } else if ((hits_vec[zoomlevel + 1]->size() <= max_hits_pairs) &&
547  ((zoomlevel + 1) >= min_zoom)) {
548  findHelicesByPairsBegin(min_hits, max_hits, tracks, maxtracks,
549  zoomlevel + 1);
550  } else {
551  findHelices(min_hits, max_hits, tracks, maxtracks, zoomlevel + 1);
552  }
553  if (maxtracks != 0) {
554  double phi_proportion =
555  (zoomranges[zoomlevel].max_phi -
556  zoomranges[zoomlevel].min_phi) /
557  ((zoomranges[0].max_phi - zoomranges[0].min_phi));
558  double d_proportion =
559  (zoomranges[zoomlevel].max_d - zoomranges[zoomlevel].min_d) /
560  ((zoomranges[0].max_d - zoomranges[0].min_d));
561  double k_proportion =
562  (zoomranges[zoomlevel].max_k - zoomranges[zoomlevel].min_k) /
563  ((zoomranges[0].max_k - zoomranges[0].min_k));
564  unsigned int expected_tracks =
565  (unsigned int)(fabs(((double)maxtracks) * phi_proportion *
566  d_proportion * k_proportion)) +
567  1;
568 
569  if ((tracks.size() - tracks_at_start) > expected_tracks) {
570  return;
571  }
572  }
573  }
574  if (count == n_entries) {
575  break;
576  }
577  hits_vec[zoomlevel + 1]->clear();
578  setRange((*(bins_vec[zoomlevel]))[count], zoomranges[zoomlevel],
579  zoomranges[zoomlevel + 1], n_phi, n_d, n_k, n_dzdl, n_z0);
580  }
581  }
582  } else {
583  if (clusters_vec[zoomlevel]->size() == 0) {
584  return;
585  }
586  // for each cluster, eiter perform the Hough again or break into
587  // user-defined routine
588  for (unsigned int i = 0, size = num_clusters[zoomlevel]; i < size; ++i) {
589  hits_vec[zoomlevel + 1]->clear();
590  vector<unsigned int>::iterator index_iter;
591  for (index_iter = (*(clusters_vec[zoomlevel]))[i].hit_indexes.begin();
592  index_iter != (*(clusters_vec[zoomlevel]))[i].hit_indexes.end();
593  ++index_iter) {
594  if (remove_hits == true) {
595  if ((*hit_used)[(*(hits_vec[zoomlevel]))[*index_iter].get_id()] ==
596  false) {
597  hits_vec[zoomlevel + 1]->push_back(
598  (*(hits_vec[zoomlevel]))[*index_iter]);
599  }
600  } else {
601  hits_vec[zoomlevel + 1]->push_back(
602  (*(hits_vec[zoomlevel]))[*index_iter]);
603  }
604  }
605 
606  setClusterRange(zoomranges[zoomlevel], zoomranges[zoomlevel + 1],
607  (*(clusters_vec[zoomlevel]))[i].range, n_phi, n_d, n_k,
608  n_dzdl, n_z0);
609  if ((breakRecursion(*(hits_vec[zoomlevel + 1]),
610  zoomranges[zoomlevel + 1]) == true)) {
611  } else if ((zoomlevel + 1) == max_zoom) {
612  findTracks(*(hits_vec[zoomlevel + 1]), tracks,
613  zoomranges[zoomlevel + 1]);
614  } else if (((zoomlevel + 1) >= min_zoom) &&
615  (hits_vec[zoomlevel + 1]->size() <= max_hits)) {
616  findTracks(*(hits_vec[zoomlevel + 1]), tracks,
617  zoomranges[zoomlevel + 1]);
618  } else if ((hits_vec[zoomlevel + 1]->size() <= max_hits_pairs) &&
619  ((zoomlevel + 1) >= min_zoom)) {
620  findHelicesByPairsBegin(min_hits, max_hits, tracks, maxtracks,
621  zoomlevel + 1);
622  } else {
623  findHelices(min_hits, max_hits, tracks, maxtracks, zoomlevel + 1);
624  }
625  }
626  }
627 }
628 
629 void HelixHough::findSeededHelices(vector<SimpleTrack3D>& seeds,
630  vector<SimpleHit3D>& hits,
631  unsigned int min_hits, unsigned int max_hits,
632  vector<SimpleTrack3D>& tracks,
633  unsigned int maxtracks) {
634  unsigned int n_layers_orig = n_layers;
635  unsigned int req_layers_orig = req_layers;
636 
637  bool smooth_back_orig = smooth_back;
638 
639  int max_layer = 0;
640  int min_layer = 999999;
641  for (unsigned int i = 0; i < hits.size(); ++i) {
642  if (hits[i].get_layer() > max_layer) {
643  max_layer = hits[i].get_layer();
644  }
645  if (hits[i].get_layer() < min_layer) {
646  min_layer = hits[i].get_layer();
647  }
648  if (hits[i].get_layer() < 0) {
649  findSeededHelices_run(seeds, hits, min_hits, max_hits, tracks, maxtracks);
650  return;
651  }
652  }
653 
654  if (min_layer == 999999) {
655  findSeededHelices_run(seeds, hits, min_hits, max_hits, tracks, maxtracks);
656  return;
657  }
658 
659  n_layers = min_layer;
660 
661  float req_prop = ((float)req_layers) / ((float)((1 + max_layer - min_layer)));
662 
663  int layers_left = (1 + max_layer - min_layer);
664  unsigned int cur_layer = min_layer;
665  vector<SimpleTrack3D> cur_seeds = seeds;
666  vector<SimpleTrack3D> cur_tracks;
667  vector<SimpleHit3D> cur_hits;
668  smooth_back = false;
669  while (layers_left > 0) {
670  unsigned int layers_iter = layers_at_a_time;
671  if ((cur_layer + layers_at_a_time - 1) > (unsigned int) max_layer) {
672  layers_iter = (max_layer - (cur_layer - 1));
673  }
674 
675  cur_tracks.clear();
676  cur_hits.clear();
677  for (unsigned int i = 0; i < hits.size(); ++i) {
678  if ((hits[i].get_layer() >= (int) cur_layer) &&
679  (hits[i].get_layer() <= (int) (cur_layer + layers_at_a_time - 1))) {
680  cur_hits.push_back(hits[i]);
681  }
682  }
683 
684  n_layers += layers_iter;
685 
686  unsigned int nreq = floor(((float)(layers_iter)) * req_prop);
687  if (nreq == 0) {
688  nreq = 1;
689  }
690  requireLayers(nreq);
691 
692  clear();
693  req_layers = nreq;
694 
695  if (layers_left <= (int) layers_at_a_time) {
696  smooth_back = smooth_back_orig;
697  }
698 
699  findSeededHelices_run(cur_seeds, cur_hits, nreq, max_hits, cur_tracks,
700  maxtracks);
701  setSeedStates(track_states);
702  cur_seeds = cur_tracks;
703  layers_left -= layers_at_a_time;
704  cur_layer += layers_at_a_time;
705  }
706  tracks = cur_tracks;
707 
708  n_layers = n_layers_orig;
709  req_layers = req_layers_orig;
710 }
711 
712 void HelixHough::findSeededHelices_run(vector<SimpleTrack3D>& seeds,
713  vector<SimpleHit3D>& hits,
714  unsigned int min_hits,
715  unsigned int max_hits,
716  vector<SimpleTrack3D>& tracks,
717  unsigned int maxtracks) {
718  vote_time = 0.;
719  xy_vote_time = 0.;
720  z_vote_time = 0.;
721 
722  initEvent(hits, min_hits);
723  initSeeding(seeds);
724 
725  hit_used->clear();
726  index_mapping.clear();
727  index_mapping.resize(hits.size(), 0);
728  for (unsigned int i = 0; i < hits.size(); i++) {
729  index_mapping[i] = hits[i].get_id();
730  hits[i].set_id(i);
731  }
732  (*hit_used).assign(hits.size(), false);
733 
734  max_tracks = maxtracks;
735  base_hits = &hits;
736  (*(hits_vec[0])) = hits;
737  (*(seeds_vec[0])) = seeds;
738  current_range = top_range;
739  zoomranges.clear();
740  for (unsigned int z = 0; z <= max_zoom; z++) {
741  zoomranges.push_back(top_range);
742  }
743  vector<SimpleTrack3D> temp_tracks;
744 
745  if (separate_by_helicity == true) {
746  helicity = true;
747  findSeededHelices(min_hits, max_hits, temp_tracks, maxtracks, 0);
748  helicity = false;
749  findSeededHelices(min_hits, max_hits, temp_tracks, maxtracks, 0);
750  } else {
751  findSeededHelices(min_hits, max_hits, temp_tracks, maxtracks, 0);
752  }
753 
754  for (unsigned int i = 0; i < hits.size(); i++) {
755  hits[i].set_id(index_mapping[i]);
756  }
757 
758  finalize(temp_tracks, tracks);
759 
760  if (print_timings == true) {
761  cout << "vote time = " << vote_time << endl;
762  cout << "xy vote time = " << xy_vote_time << endl;
763  cout << "z vote time = " << z_vote_time << endl;
764  cout << endl;
765  }
766 }
767 
768 class floatBin {
769  public:
770  floatBin(float l, float h) : low(l), high(h) {}
771  ~floatBin() {}
772  float low;
773  float high;
774  bool operator<(const floatBin& other) const { return (high < other.low); }
775 };
776 
777 // return which bin in bins that val belongs to, or -1 if it doesn't belong to
778 // any
779 static int seed_bin(vector<floatBin>& bins, float val) {
780  floatBin bin(val, val);
781  if ((bin < bins[0]) || (bins.back() < bin)) {
782  return -1;
783  }
784  pair<vector<floatBin>::iterator, vector<floatBin>::iterator> bounds;
785  bounds = equal_range(bins.begin(), bins.end(), bin);
786  return ((int)(bounds.first - bins.begin()));
787 }
788 
789 void HelixHough::setRangeFromSeed(HelixRange& range, SimpleTrack3D& seed) {
790  range.min_phi = seed.phi - 0.001;
791  range.max_phi = seed.phi + 0.001;
792  if (range.min_phi < 0.) {
793  range.min_phi = 0.;
794  }
795  if (range.max_phi > 2. * M_PI) {
796  range.max_phi = 2. * M_PI;
797  }
798  range.min_d = seed.d - 0.001;
799  range.max_d = seed.d + 0.001;
800  range.min_k = seed.kappa - 0.00001;
801  range.max_k = seed.kappa + 0.00001;
802  if (range.min_k < 0.) {
803  range.min_k = 0.;
804  }
805 
806  range.min_k = range.min_k * range.min_k;
807  range.max_k = range.max_k * range.max_k;
808 
809  range.min_dzdl = seed.dzdl - 0.05;
810  range.max_dzdl = seed.dzdl + 0.05;
811  range.min_z0 = seed.z0 - 0.01;
812  range.max_z0 = seed.z0 + 0.01;
813 }
814 
815 void HelixHough::findSeededHelices(unsigned int min_hits, unsigned int max_hits,
816  vector<SimpleTrack3D>& tracks,
817  unsigned int maxtracks,
818  unsigned int zoomlevel) {
819  unsigned int tracks_at_start = tracks.size();
820 
821  if ((maxtracks != 0) && (tracks.size() >= max_tracks)) {
822  return;
823  }
824 
825  timeval t1, t2;
826  double time1 = 0.;
827  double time2 = 0.;
828  gettimeofday(&t1, NULL);
829  vote(zoomlevel);
830  gettimeofday(&t2, NULL);
831  time1 = ((double)(t1.tv_sec) + (double)(t1.tv_usec) / 1000000.);
832  time2 = ((double)(t2.tv_sec) + (double)(t2.tv_usec) / 1000000.);
833  vote_time += (time2 - time1);
834 
835  unsigned int n_phi = n_phi_bins[zoomlevel];
836  unsigned int n_d = n_d_bins[zoomlevel];
837  unsigned int n_k = n_k_bins[zoomlevel];
838  unsigned int n_dzdl = n_dzdl_bins[zoomlevel];
839  unsigned int n_z0 = n_z0_bins[zoomlevel];
840 
841  float low = 0.;
842  float high = 0.;
843  float delta = 0.;
844 
845  vector<floatBin> phibins;
846  low = zoomranges[zoomlevel].min_phi;
847  delta = (zoomranges[zoomlevel].max_phi - zoomranges[zoomlevel].min_phi) /
848  ((float)(n_phi));
849  high = low + delta;
850  for (unsigned int i = 0; i < n_phi; ++i) {
851  phibins.push_back(floatBin(low, high));
852  low += delta;
853  high += delta;
854  }
855 
856  vector<floatBin> dbins;
857  low = zoomranges[zoomlevel].min_d;
858  delta = (zoomranges[zoomlevel].max_d - zoomranges[zoomlevel].min_d) /
859  ((float)(n_d));
860  high = low + delta;
861  for (unsigned int i = 0; i < n_d; ++i) {
862  dbins.push_back(floatBin(low, high));
863  low += delta;
864  high += delta;
865  }
866 
867  vector<floatBin> kbins;
868  low = zoomranges[zoomlevel].min_k;
869  delta = (zoomranges[zoomlevel].max_k - zoomranges[zoomlevel].min_k) /
870  ((float)(n_k));
871  high = low + delta;
872  for (unsigned int i = 0; i < n_k; ++i) {
873  kbins.push_back(floatBin(low, high));
874  low += delta;
875  high += delta;
876  }
877 
878  vector<floatBin> z0bins;
879  low = zoomranges[zoomlevel].min_z0;
880  delta = (zoomranges[zoomlevel].max_z0 - zoomranges[zoomlevel].min_z0) /
881  ((float)(n_z0));
882  high = low + delta;
883  for (unsigned int i = 0; i < n_z0; ++i) {
884  z0bins.push_back(floatBin(low, high));
885  low += delta;
886  high += delta;
887  }
888 
889  vector<floatBin> dzdlbins;
890  low = zoomranges[zoomlevel].min_dzdl;
891  delta = (zoomranges[zoomlevel].max_dzdl - zoomranges[zoomlevel].min_dzdl) /
892  ((float)(n_dzdl));
893  high = low + delta;
894  for (unsigned int i = 0; i < n_dzdl; ++i) {
895  dzdlbins.push_back(floatBin(low, high));
896  low += delta;
897  high += delta;
898  }
899 
900  // voting for the seeds
901  for (unsigned int i = 0; i < seeds_vec[zoomlevel]->size(); ++i) {
902  float d = (*(seeds_vec[zoomlevel]))[i].d;
903  float phi = (*(seeds_vec[zoomlevel]))[i].phi;
904  float kappa = pow((*(seeds_vec[zoomlevel]))[i].kappa, 1.);
905  float z0 = (*(seeds_vec[zoomlevel]))[i].z0;
906  float dzdl = (*(seeds_vec[zoomlevel]))[i].dzdl;
907 
908  int tempbin = 0;
909 
910  unsigned int phi_bin = 0;
911  tempbin = seed_bin(phibins, phi);
912  if (tempbin < 0) {
913  continue;
914  } else {
915  phi_bin = (unsigned int)(tempbin);
916  }
917 
918  unsigned int d_bin = 0;
919  tempbin = seed_bin(dbins, d);
920  if (tempbin < 0) {
921  continue;
922  } else {
923  d_bin = (unsigned int)(tempbin);
924  }
925 
926  unsigned int kappa_bin = 0;
927  tempbin = seed_bin(kbins, kappa);
928  if (tempbin < 0) {
929  continue;
930  } else {
931  kappa_bin = (unsigned int)(tempbin);
932  }
933 
934  unsigned int z0_bin = 0;
935  tempbin = seed_bin(z0bins, z0);
936  if (tempbin < 0) {
937  continue;
938  } else {
939  z0_bin = (unsigned int)(tempbin);
940  }
941 
942  unsigned int dzdl_bin = 0;
943  tempbin = seed_bin(dzdlbins, dzdl);
944  if (tempbin < 0) {
945  continue;
946  } else {
947  dzdl_bin = (unsigned int)(tempbin);
948  }
949 
950  unsigned int bin = BinEntryPair5D::linearBin(
951  n_d, n_k, n_dzdl, n_z0, phi_bin, d_bin, kappa_bin, dzdl_bin, z0_bin);
952 
953  bins_vec[zoomlevel]->push_back(BinEntryPair5D(bin, i, true));
954  }
955  sort(bins_vec[zoomlevel]->begin(), bins_vec[zoomlevel]->end());
956 
957  unsigned int n_entries = bins_vec[zoomlevel]->size();
958  if (n_entries == 0) {
959  return;
960  }
961 
962  unsigned int count = 0;
963  hits_vec[zoomlevel + 1]->clear();
964  seeds_vec[zoomlevel + 1]->clear();
965  setRange((*(bins_vec[zoomlevel]))[count], zoomranges[zoomlevel],
966  zoomranges[zoomlevel + 1], n_phi, n_d, n_k, n_dzdl, n_z0);
967  // scan over the bins in 5-D hough space
968  while (count < n_entries) {
969  if ((*(bins_vec[zoomlevel]))[count].is_seed == false) {
970  if (remove_hits == true) {
971  if ((*hit_used)[(*(
972  hits_vec[zoomlevel]))[(*(bins_vec[zoomlevel]))[count].entry]
973  .get_id()] == false) {
974  hits_vec[zoomlevel + 1]->push_back(
975  (*(hits_vec[zoomlevel]))[(*(bins_vec[zoomlevel]))[count].entry]);
976  }
977  } else {
978  hits_vec[zoomlevel + 1]->push_back(
979  (*(hits_vec[zoomlevel]))[(*(bins_vec[zoomlevel]))[count].entry]);
980  }
981  } else {
982  seeds_vec[zoomlevel + 1]->push_back(
983  (*(seeds_vec[zoomlevel]))[(*(bins_vec[zoomlevel]))[count].entry]);
984  }
985 
986  count++;
987 
988  // we have collected all hits from this bin. now zoom again or find tracks
989  // with user routine
990  if ((count == n_entries) || ((*(bins_vec[zoomlevel]))[count].bin !=
991  (*(bins_vec[zoomlevel]))[count - 1].bin)) {
992  if ((hits_vec[zoomlevel + 1]->size() >= min_hits) &&
993  (seeds_vec[zoomlevel + 1]->size() != 0)) {
994  if (breakRecursion(*(hits_vec[zoomlevel + 1]),
995  zoomranges[zoomlevel + 1]) == true) {
996  findSeededTracks(*(seeds_vec[zoomlevel + 1]),
997  *(hits_vec[zoomlevel + 1]), tracks,
998  zoomranges[zoomlevel + 1]);
999  } else if (((zoomlevel + 1) >= min_zoom) &&
1000  ((hits_vec[zoomlevel + 1]->size() <= max_hits) ||
1001  (zoomlevel + 1) == max_zoom)) {
1002  findSeededTracks(*(seeds_vec[zoomlevel + 1]),
1003  *(hits_vec[zoomlevel + 1]), tracks,
1004  zoomranges[zoomlevel + 1]);
1005  } else {
1006  findSeededHelices(min_hits, max_hits, tracks, maxtracks,
1007  zoomlevel + 1);
1008  }
1009  if (maxtracks != 0) {
1010  double phi_proportion =
1011  (zoomranges[zoomlevel].max_phi - zoomranges[zoomlevel].min_phi) /
1012  ((zoomranges[0].max_phi - zoomranges[0].min_phi));
1013  double d_proportion =
1014  (zoomranges[zoomlevel].max_d - zoomranges[zoomlevel].min_d) /
1015  ((zoomranges[0].max_d - zoomranges[0].min_d));
1016  double k_proportion =
1017  (zoomranges[zoomlevel].max_k - zoomranges[zoomlevel].min_k) /
1018  ((zoomranges[0].max_k - zoomranges[0].min_k));
1019  unsigned int expected_tracks =
1020  (unsigned int)(fabs(((double)maxtracks) * phi_proportion *
1021  d_proportion * k_proportion)) +
1022  1;
1023 
1024  if ((tracks.size() - tracks_at_start) > expected_tracks) {
1025  return;
1026  }
1027  }
1028  }
1029  if (count == n_entries) {
1030  break;
1031  }
1032  hits_vec[zoomlevel + 1]->clear();
1033  seeds_vec[zoomlevel + 1]->clear();
1034  setRange((*(bins_vec[zoomlevel]))[count], zoomranges[zoomlevel],
1035  zoomranges[zoomlevel + 1], n_phi, n_d, n_k, n_dzdl, n_z0);
1036  }
1037  }
1038 }