ECCE @ EIC Software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
sPHENIXTrackerTpc.h
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file sPHENIXTrackerTpc.h
1 #ifndef SPHENIX_SPHENIXTRACKERTPC_H
2 #define SPHENIX_SPHENIXTRACKERTPC_H
3 
4 #include "HelixHough.h"
5 #include "HelixKalmanState.h"
6 #include "HelixRange.h"
7 #include "Seamstress.h"
8 #include "SimpleHit3D.h"
9 
10 #include <algorithm>
11 #include <cmath>
12 #include <map>
13 #include <memory>
14 #include <set>
15 #include <vector>
16 
17 class CylinderKalman;
18 class HelixResolution;
19 class SimpleTrack3D;
20 
21 namespace SeamStress { template <class TClass> class Pincushion; }
22 
23 class AngleIndexPair {
24  public:
25  AngleIndexPair(float ang, unsigned int idx) : angle(ang), index(idx) {
26  float twopi = 2. * M_PI;
27  int a = (int)(angle / twopi);
28  angle -= a * twopi;
29  while (angle < 0.) {
30  angle += twopi;
31  }
32  while (angle >= twopi) {
33  angle -= twopi;
34  }
35  }
37 
38  bool operator<(const AngleIndexPair& other) const {
39  return angle < other.angle;
40  }
41 
42  static float absDiff(float angle1, float angle2) {
43  float diff = (angle1 - angle2);
44  while (diff > M_PI) {
45  diff -= 2. * M_PI;
46  }
47  while (diff < -M_PI) {
48  diff += 2. * M_PI;
49  }
50  diff = fabs(diff);
51  return diff;
52  }
53 
54  float angle;
55  unsigned int index;
56 };
57 
58 class AngleIndexList {
59  public:
62 
63  void addPair(AngleIndexPair& angind) {
64  sorted = false;
65  vec.push_back(angind);
66  }
67 
68  void getRangeListSimple(float angle, float error,
69  std::vector<AngleIndexPair*>& result) {
70  result.clear();
71 
72  for (unsigned int i = 0; i < vec.size(); i++) {
73  if (AngleIndexPair::absDiff(angle, vec[i].angle) <= error) {
74  result.push_back(&(vec[i]));
75  }
76  }
77  }
78 
79  void getRangeList(float angle, float error,
80  std::vector<AngleIndexPair*>& result) {
81  float twopi = 2. * M_PI;
82  int a = (int)(angle / twopi);
83  angle -= a * twopi;
84  while (angle < 0.) {
85  angle += twopi;
86  }
87  while (angle >= twopi) {
88  angle -= twopi;
89  }
90 
91  if (vec.size() <= 4) {
92  return getRangeListSimple(angle, error, result);
93  }
94 
95  result.clear();
96 
97  unsigned int closest = findClosest(angle);
98  // first, traverse upward
99  unsigned int current = closest;
100  unsigned int lowest = 0;
101  unsigned int highest = vec.size() - 1;
102  while (true) {
103  if (AngleIndexPair::absDiff(angle, vec[current].angle) <= error) {
104  result.push_back(&(vec[current]));
105  current = (current + 1) % (vec.size());
106  if (current == closest) {
107  return;
108  }
109  } else {
110  break;
111  }
112  }
113 
114  // now, traverse downward
115  if (current <= closest) {
116  lowest = current;
117  } else {
118  highest = current;
119  }
120  current = ((closest + vec.size()) - 1) % (vec.size());
121  while (true) {
122  if ((current == lowest) || (current == highest)) {
123  break;
124  }
125  if (AngleIndexPair::absDiff(angle, vec[current].angle) <= error) {
126  result.push_back(&(vec[current]));
127  current = ((current + vec.size()) - 1) % (vec.size());
128  } else {
129  break;
130  }
131  }
132  }
133 
134  unsigned int findClosestSimple(float angle, unsigned int lower,
135  unsigned int upper) {
136  unsigned int closest = lower;
137  float diff = AngleIndexPair::absDiff(vec[closest].angle, angle);
138  for (unsigned int i = (lower + 1); i <= upper; i++) {
139  float tempdiff = AngleIndexPair::absDiff(vec[i].angle, angle);
140  if (tempdiff < diff) {
141  closest = i;
142  diff = tempdiff;
143  }
144  }
145 
146  return closest;
147  }
148 
149  unsigned int findClosest(float angle) {
150  if (vec.size() <= 4) {
151  return findClosestSimple(angle, 0, vec.size() - 1);
152  }
153 
154  if (sorted == false) {
155  std::sort(vec.begin(), vec.end());
156  sorted = true;
157  }
158 
159  unsigned int lower = 0;
160  unsigned int upper = vec.size() - 1;
161  unsigned int middle = vec.size() / 2;
162  while (true) {
163  if ((upper - lower) <= 4) {
164  return findClosestSimple(angle, lower, upper);
165  }
166 
167  if (angle <= vec[middle].angle) {
168  upper = middle;
169  middle = (lower + upper) / 2;
170  } else {
171  lower = middle;
172  middle = (lower + upper) / 2;
173  }
174  }
175  }
176 
177  std::vector<AngleIndexPair> vec;
178  bool sorted;
179 };
180 
181 class TrackSegment {
182  public:
184  : chi2(0.), ux(0.), uy(0.), kappa(0.), dkappa(0.), seed(0), n_hits(0) {}
186 
187  float chi2;
188  float ux, uy;
189  float kappa;
190  float dkappa;
191  std::vector<unsigned int> hits;
192  unsigned int seed;
193  unsigned int bin;
194  unsigned int n_hits;
195 };
196 
198  public:
199  sPHENIXTrackerTpc(unsigned int n_phi, unsigned int n_d, unsigned int n_k,
200  unsigned int n_dzdl, unsigned int n_z0,
201  HelixResolution& min_resolution,
202  HelixResolution& max_resolution, HelixRange& range,
203  std::vector<float>& material, std::vector<float>& radius,
204  float Bfield);
205  sPHENIXTrackerTpc(std::vector<std::vector<unsigned int> >& zoom_profile,
206  unsigned int minzoom, HelixRange& range,
207  std::vector<float>& material, std::vector<float>& radius,
208  float Bfield, bool parallel = false,
209  unsigned int num_threads = 1);
210  virtual ~sPHENIXTrackerTpc();
211 
212  void finalize(std::vector<SimpleTrack3D>& input,
213  std::vector<SimpleTrack3D>& output);
214 
215  void findTracks(std::vector<SimpleHit3D>& hits,
216  std::vector<SimpleTrack3D>& tracks,
217  const HelixRange& range);
218  void findTracksBySegments(std::vector<SimpleHit3D>& hits,
219  std::vector<SimpleTrack3D>& tracks,
220  const HelixRange& range);
221  void initEvent(std::vector<SimpleHit3D>& hits, unsigned int /*min_hits*/) {
222 
223  int min_layer = 999999;
224  int max_layer = 0;
225  for (unsigned int i = 0; i < hits.size(); ++i) {
226  if (hits[i].get_layer() < min_layer) {
227  min_layer = hits[i].get_layer();
228  }
229  if (hits[i].get_layer() > max_layer) {
230  max_layer = hits[i].get_layer();
231  }
232  }
233  setSeedLayer(min_layer);
234  setNLayers(max_layer + 1);
235 
236 
237  nfits = 0;
238  combos.clear();
239  fail_combos.clear();
240  pass_combos.clear();
241  seeding = false;
242  // required_layers = min_hits;
244  findtracksiter = 0;
245  CAtime = 0.;
246  KALtime = 0.;
247  findtracks_bin = 0;
248  }
249 
250  void initSeeding(std::vector<SimpleTrack3D>& seeds) {
251  seeding = true;
252  seed_used.clear();
253  seed_used.assign(seeds.size(), false);
254  AngleIndexList templist;
255  angle_list.clear();
256  angle_list.assign(n_layers, templist);
257  }
258 
259  void clear() {
260  track_states.clear();
261  isolation_variable.clear();
262  }
263 
264  bool breakRecursion(const std::vector<SimpleHit3D>& hits,
265  const HelixRange& range);
266  float phiError(SimpleHit3D& hit, float min_k, float max_k, float min_d,
267  float max_d, float min_z0, float max_z0, float min_dzdl,
268  float max_dzdl, bool pairvoting = false);
269  float dzdlError(SimpleHit3D& hit, float min_k, float max_k, float min_d,
270  float max_d, float min_z0, float max_z0, float min_dzdl,
271  float max_dzdl, bool pairvoting = false);
272 
273  static float fitTrack(SimpleTrack3D& track,
274  float scale = 1.0);
275  static float fitTrack(SimpleTrack3D& track,
276  std::vector<float>& chi2_hit,
277  float scale = 1.0);
278 
279  void setVerbosity(int v) { verbosity = v; }
280 
281  void setCutOnDca(bool dcut) { cut_on_dca = dcut; }
282  void setDcaCut(float dcut) { dca_cut = dcut; }
283  void setVertex(float vx, float vy, float vz) {
284  vertex_x = vx;
285  vertex_y = vy;
286  vertex_z = vz;
287  }
288 
289  void setRejectGhosts(bool rg) { reject_ghosts = rg; }
290 
291  std::vector<float>& getIsolation() { return isolation_variable; }
292 
293  static void calculateKappaTangents(
294  float* x1_a, float* y1_a, float* z1_a, float* x2_a, float* y2_a,
295  float* z2_a, float* x3_a, float* y3_a, float* z3_a, float* dx1_a,
296  float* dy1_a, float* dz1_a, float* dx2_a, float* dy2_a, float* dz2_a,
297  float* dx3_a, float* dy3_a, float* dz3_a, float* kappa_a, float* dkappa_a,
298  float* ux_mid_a, float* uy_mid_a, float* ux_end_a, float* uy_end_a,
299  float* dzdl_1_a, float* dzdl_2_a, float* ddzdl_1_a, float* ddzdl_2_a);
301  float* x1_a, float* y1_a, float* z1_a, float* x2_a, float* y2_a,
302  float* z2_a, float* x3_a, float* y3_a, float* z3_a, float* dx1_a,
303  float* dy1_a, float* dz1_a, float* dx2_a, float* dy2_a, float* dz2_a,
304  float* dx3_a, float* dy3_a, float* dz3_a, float* kappa_a, float* dkappa_a,
305  float* ux_mid_a, float* uy_mid_a, float* ux_end_a, float* uy_end_a,
306  float* dzdl_1_a, float* dzdl_2_a, float* ddzdl_1_a, float* ddzdl_2_a,
307  float sinang_cut, float cosang_diff_inv, float* cur_kappa_a,
308  float* cur_dkappa_a, float* cur_ux_a, float* cur_uy_a, float* cur_chi2_a,
309  float* chi2_a);
310 
311  void pairRejection(std::vector<SimpleTrack3D>& input,
312  std::vector<SimpleTrack3D>& output,
313  std::vector<bool>& usetrack,
314  std::vector<float>& next_best_chi2);
315  void tripletRejection(std::vector<SimpleTrack3D>& input,
316  std::vector<SimpleTrack3D>& output,
317  std::vector<bool>& usetrack,
318  std::vector<float>& next_best_chi2);
319 
320  bool seedWasUsed(unsigned int seed_index) { return seed_used.at(seed_index); }
321 
322  void setSeparateByHelicity(bool sbh) {
323  separate_by_helicity = sbh;
324  if (is_parallel == true) {
325  for (unsigned int i = 0; i < thread_trackers.size(); ++i) {
326  thread_trackers[i]->setSeparateByHelicity(sbh);
327  }
328  }
329  }
330 
331  void requireLayers(unsigned int nl) {
332  check_layers = true;
333  req_layers = nl;
334  if (is_parallel == true) {
335  for (unsigned int i = 0; i < thread_trackers.size(); ++i) {
336  thread_trackers[i]->requireLayers(nl);
337  }
338  }
339  }
340 
341  virtual void setMaxHitsPairs(unsigned int mhp) {
342  max_hits_pairs = mhp;
343  if (is_parallel == true) {
344  for (unsigned int i = 0; i < thread_trackers.size(); ++i) {
345  thread_trackers[i]->setMaxHitsPairs(mhp);
346  }
347  }
348  }
349 
350  void setBinScale(float b_scl) {
351  bin_scale = b_scl;
352  if (is_parallel == true) {
353  for (unsigned int i = 0; i < thread_trackers.size(); ++i) {
354  thread_trackers[i]->setBinScale(b_scl);
355  }
356  }
357  }
358  void setZBinScale(float b_scl) {
359  z_bin_scale = b_scl;
360  if (is_parallel == true) {
361  for (unsigned int i = 0; i < thread_trackers.size(); ++i) {
362  thread_trackers[i]->setZBinScale(b_scl);
363  }
364  }
365  }
366  void setRemoveHits(bool rh) {
367  remove_hits = rh;
368  if (is_parallel == true) {
369  for (unsigned int i = 0; i < thread_trackers.size(); ++i) {
370  thread_trackers[i]->setRemoveHits(rh);
371  }
372  }
373  }
374 
375 
376  void setStartLayer(int start) {
377  layer_start = start;
378  }
379 
380  void setEndLayer(int end) {
381  layer_end = end;
382  }
383 
384  void setSeedLayer(int sl) {
385  seed_layer = sl;
386  if (is_parallel == true) {
387  for (unsigned int i = 0; i < thread_trackers.size(); ++i) {
388  thread_trackers[i]->setSeedLayer(sl);
389  }
390  }
391  }
392 
393  void setNLayers(unsigned int n) {
394  n_layers = n;
395  std::vector<SimpleHit3D> one_layer;
396  layer_sorted.clear();
397  layer_sorted.assign(n_layers, one_layer);
398  temp_comb.assign(n_layers, 0);
399  if (is_parallel == true) {
400  for (unsigned int i = 0; i < thread_trackers.size(); ++i) {
401  thread_trackers[i]->setNLayers(n);
402  }
403  }
404  }
405 
406  void setFastChi2Cut(float par0, float par1, float max) {
407  fast_chi2_cut_par0 = par0;
408  fast_chi2_cut_par1 = par1;
410  if (is_parallel == true) {
411  for (unsigned int i = 0; i < thread_trackers.size(); ++i) {
412  thread_trackers[i]->setFastChi2Cut(par0, par1, max);
413  }
414  }
415  }
416  void setChi2Cut(float c) {
417  chi2_cut = c;
418  if (is_parallel == true) {
419  for (unsigned int i = 0; i < thread_trackers.size(); ++i) {
420  thread_trackers[i]->setChi2Cut(c);
421  }
422  }
423  }
424  void setChi2RemovalCut(float c) {
426  if (is_parallel == true) {
427  for (unsigned int i = 0; i < thread_trackers.size(); ++i) {
428  thread_trackers[i]->setChi2RemovalCut(c);
429  }
430  }
431  }
432 
433  void setNRemovalHits(unsigned int n) {
434  n_removal_hits = n;
435  if (is_parallel == true) {
436  for (unsigned int i = 0; i < thread_trackers.size(); ++i) {
437  thread_trackers[i]->setNRemovalHits(n);
438  }
439  }
440  }
441 
442  void setCosAngleCut(float cut) {
443  cosang_cut = cut;
444  if (is_parallel == true) {
445  for (unsigned int i = 0; i < thread_trackers.size(); ++i) {
446  thread_trackers[i]->setCosAngleCut(cut);
447  }
448  }
449  }
450 
451  void setCellularAutomatonChi2Cut(float cut) {
452  ca_chi2_cut = cut;
453  if (is_parallel == true) {
454  for (unsigned int i = 0; i < thread_trackers.size(); ++i) {
455  thread_trackers[i]->setCellularAutomatonChi2Cut(cut);
456  }
457  }
458  }
459 
460  void setThread() { is_thread = true; }
461 
462  void initSplitting(std::vector<SimpleHit3D>& hits, unsigned int min_hits,
463  unsigned int max_hits);
464 
465  void setHitErrorScale(unsigned int layer, float scale) {
466  if (layer >= hit_error_scale.size()) {
467  hit_error_scale.resize(layer + 1, 1.);
468  }
470  }
471 
472  void setRequirePixels(bool rp){require_pixels = rp;}
473 
474 
475  private:
476  float kappaToPt(float kappa);
477  float ptToKappa(float pt);
478  void findTracksByCombinatorialKalman(std::vector<SimpleHit3D>& hits,
479  std::vector<SimpleTrack3D>& tracks,
480  const HelixRange& range);
481 
485  float chi2_cut;
487  unsigned int n_removal_hits;
488  std::set<std::vector<unsigned int> > combos;
489  std::map<std::vector<unsigned int>, HelixKalmanState> pass_combos;
490  std::set<std::vector<unsigned int> > fail_combos;
491  std::vector<float> detector_scatter;
492  std::vector<float> integrated_scatter;
493  std::vector<float> detector_radii;
494  bool seeding;
497  float dca_cut;
499  unsigned int required_layers;
501  unsigned int nfits;
502 
504 
505  std::vector<float> isolation_variable;
506 
507  unsigned int findtracksiter;
508 
509  float prev_max_k;
511  float prev_p_inv;
512 
514 
515  std::vector<TrackSegment> segments1;
516  std::vector<TrackSegment> segments2;
517 
518  std::vector<std::vector<SimpleHit3D> > layer_sorted;
519  std::vector<unsigned int> temp_comb;
520 
522 
523  std::vector<bool> seed_used;
524 
525  unsigned int nthreads;
526  std::vector<SeamStress::Seamstress*> *vssp;
527  std::vector<SeamStress::Seamstress> vss;
529  std::vector<sPHENIXTrackerTpc*> thread_trackers;
530  std::vector<std::vector<SimpleTrack3D> > thread_tracks;
531  std::vector<HelixRange> thread_ranges;
532  std::vector<std::vector<SimpleHit3D> > thread_hits;
533  std::vector<std::vector<SimpleHit3D> > split_input_hits;
534  std::vector<std::vector<std::vector<SimpleHit3D> >* > split_output_hits;
535  std::vector<std::vector<HelixRange>* > split_ranges;
537  unsigned int thread_min_hits;
538  unsigned int thread_max_hits;
539  bool is_thread;
540 
541  std::vector<AngleIndexList> angle_list;
542 
543  double CAtime;
544  double KALtime;
545 
546  std::vector<std::vector<SimpleHit3D> > layer_sorted_1[4];
547  unsigned int findtracks_bin;
548 
549  float ca_chi2_cut;
550  float cosang_cut;
551 
552  std::vector<float> hit_error_scale;
553 
555 
556 };
557 
558 
559 #endif