ECCE @ EIC Software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
sPHENIXTracker.cpp
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file sPHENIXTracker.cpp
1 #include "sPHENIXTracker.h"
2 #include "CylinderKalman.h"
3 #include "Pincushion.h"
4 #include "SimpleTrack3D.h"
5 #include "vector_math_inline.h"
6 
7 #include <Eigen/Core>
8 #include <Eigen/Dense>
9 #include <Eigen/LU>
10 
11 #include <algorithm>
12 #include <cfloat>
13 #include <cmath>
14 #include <cstddef>
15 #include <sys/time.h>
16 #include <iostream>
17 #include <utility>
18 #include <xmmintrin.h>
19 
20 class HelixResolution;
21 
22 using namespace std;
23 using namespace Eigen;
24 using namespace SeamStress;
25 
26 class hit_triplet {
27  public:
28  hit_triplet(unsigned int h1, unsigned int h2, unsigned int h3, unsigned int t,
29  float c)
30  : hit1(h1), hit2(h2), hit3(h3), track(t), chi2(c) {}
32 
33  bool operator<(const hit_triplet& other) const {
34  return (hit1 < other.hit1) ||
35  ((hit2 < other.hit2) && (hit1 == other.hit1)) ||
36  ((hit3 < other.hit3) && (hit1 == other.hit1) &&
37  (hit2 == other.hit2));
38  }
39 
40  bool operator==(const hit_triplet& other) const {
41  return ((hit1 == other.hit1) && (hit2 == other.hit2) &&
42  (hit3 == other.hit3));
43  }
44 
45  unsigned int hit1, hit2, hit3, track;
46  float chi2;
47 };
48 
49 class hitTriplet {
50  public:
51  hitTriplet(unsigned int h1, unsigned int h2, unsigned int h3, unsigned int t,
52  float c)
53  : hit1(h1), hit2(h2), hit3(h3), track(t), chi2(c) {}
55 
56  bool operator<(const hitTriplet& other) const {
57  return (hit1 < other.hit1) ||
58  ((hit2 < other.hit2) && (hit1 == other.hit1)) ||
59  ((hit3 < other.hit3) && (hit1 == other.hit1) &&
60  (hit2 == other.hit2));
61  }
62 
63  bool operator==(const hitTriplet& other) const {
64  return ((hit1 == other.hit1) && (hit2 == other.hit2) &&
65  (hit3 == other.hit3));
66  }
67 
68  unsigned int hit1, hit2, hit3, track;
69  float chi2;
70 };
71 
72 void sPHENIXTracker::tripletRejection(vector<SimpleTrack3D>& input,
73  vector<SimpleTrack3D>& /*output*/,
74  vector<bool>& usetrack,
75  vector<float>& next_best_chi2) {
76  vector<hitTriplet> trips;
77  for (unsigned int i = 0; i < input.size(); ++i) {
78  for (unsigned int h1 = 0; h1 < input[i].hits.size(); ++h1) {
79  for (unsigned int h2 = (h1 + 1); h2 < input[i].hits.size(); ++h2) {
80  for (unsigned int h3 = (h2 + 1); h3 < input[i].hits.size(); ++h3) {
81  if (cut_on_dca == false) {
82  trips.push_back(hitTriplet(
83  input[i].hits[h1].get_id(), input[i].hits[h2].get_id(),
84  input[i].hits[h3].get_id(), i, track_states[i].chi2));
85  }
86  }
87  }
88  }
89  }
90  if (trips.size() == 0) {
91  return;
92  }
93  sort(trips.begin(), trips.end());
94  unsigned int pos = 0;
95 // unsigned int cur_h1 = trips[pos].hit1;
96 // unsigned int cur_h2 = trips[pos].hit2;
97  while (pos < trips.size()) {
98  unsigned int next_pos = pos + 1;
99  if (next_pos >= trips.size()) {
100  break;
101  }
102  while (trips[pos] == trips[next_pos]) {
103  next_pos += 1;
104  if (next_pos >= trips.size()) {
105  break;
106  }
107  }
108  if ((next_pos - pos) > 1) {
109  float best_chi2 = trips[pos].chi2;
110  float next_chi2 = trips[pos + 1].chi2;
111  unsigned int best_pos = pos;
112  for (unsigned int i = (pos + 1); i < next_pos; ++i) {
113  if (input[trips[i].track].hits.size() <
114  input[trips[best_pos].track].hits.size()) {
115  continue;
116  } else if ((input[trips[i].track].hits.size() >
117  input[trips[best_pos].track].hits.size()) ||
118  (input[trips[i].track].hits.back().get_layer() >
119  input[trips[best_pos].track].hits.back().get_layer())) {
120  next_chi2 = best_chi2;
121  best_chi2 = trips[i].chi2;
122  best_pos = i;
123  continue;
124  }
125  if ((trips[i].chi2 < best_chi2) ||
126  (usetrack[trips[best_pos].track] == false)) {
127  next_chi2 = best_chi2;
128  best_chi2 = trips[i].chi2;
129  best_pos = i;
130  } else if (trips[i].chi2 < next_chi2) {
131  next_chi2 = trips[i].chi2;
132  }
133  }
134  for (unsigned int i = pos; i < next_pos; ++i) {
135  if (i != best_pos) {
136  usetrack[trips[i].track] = false;
137  } else {
138  next_best_chi2[trips[i].track] = next_chi2;
139  }
140  }
141  }
142  pos = next_pos;
143  // cur_h1 = trips[pos].hit1;
144  // cur_h2 = trips[pos].hit2;
145  }
146 }
147 
148 sPHENIXTracker::sPHENIXTracker(unsigned int n_phi, unsigned int n_d,
149  unsigned int n_k, unsigned int n_dzdl,
150  unsigned int n_z0,
151  HelixResolution& min_resolution,
152  HelixResolution& max_resolution,
153  HelixRange& range, vector<float>& material,
154  vector<float>& radius, float Bfield)
155  : HelixHough(n_phi, n_d, n_k, n_dzdl, n_z0, min_resolution, max_resolution,
156  range),
157  fast_chi2_cut_par0(12.),
158  fast_chi2_cut_par1(0.),
159  fast_chi2_cut_max(FLT_MAX),
160  chi2_cut(3.),
161  chi2_removal_cut(1.),
162  n_removal_hits(0),
163  seeding(false),
164  verbosity(0),
165  cut_on_dca(false),
166  dca_cut(0.01),
167  vertex_x(0.),
168  vertex_y(0.),
169  vertex_z(0.),
170  required_layers(0),
171  reject_ghosts(false),
172  nfits(0),
173  findtracksiter(0),
174  prev_max_k(0.),
175  prev_max_dzdl(0.),
176  prev_p_inv(0.),
177  seed_layer(0),
178  ca_chi2_cut(2.0),
179  cosang_cut(0.985) {
180  vector<float> detector_material;
181 
182  for (unsigned int i = 0; i < radius.size(); ++i) {
183  detector_radii.push_back(radius[i]);
184  }
185  for (unsigned int i = 0; i < material.size(); ++i) {
186  detector_scatter.push_back(1.41421356237309515 * 0.0136 *
187  sqrt(3. * material[i]));
188  detector_material.push_back(3. * material[i]);
189  }
190 
191  detector_B_field = Bfield;
192 
193  integrated_scatter.assign(detector_scatter.size(), 0.);
194  float total_scatter_2 = 0.;
195  for (unsigned int l = 0; l < detector_scatter.size(); ++l) {
196  total_scatter_2 += detector_scatter[l] * detector_scatter[l];
197  integrated_scatter[l] = sqrt(total_scatter_2);
198  }
199 
200  kalman =
201  new CylinderKalman(detector_radii, detector_material, detector_B_field);
202 
203  vector<SimpleHit3D> one_layer;
204  layer_sorted.assign(n_layers, one_layer);
205  for (unsigned int i = 0; i < 4; ++i) {
206  layer_sorted_1[i].assign(n_layers, one_layer);
207  }
208  temp_comb.assign(n_layers, 0);
209 }
210 
211 sPHENIXTracker::sPHENIXTracker(vector<vector<unsigned int> >& zoom_profile,
212  unsigned int minzoom, HelixRange& range,
213  vector<float>& material, vector<float>& radius,
214  float Bfield, bool parallel,
215  unsigned int num_threads)
216  : HelixHough(zoom_profile, minzoom, range),
217  fast_chi2_cut_par0(12.),
218  fast_chi2_cut_par1(0.),
219  fast_chi2_cut_max(FLT_MAX),
220  chi2_cut(3.),
221  chi2_removal_cut(1.),
222  n_removal_hits(0),
223  seeding(false),
224  verbosity(0),
225  cut_on_dca(false),
226  dca_cut(0.01),
227  vertex_x(0.),
228  vertex_y(0.),
229  vertex_z(0.),
230  required_layers(0),
231  reject_ghosts(false),
232  nfits(0),
233  findtracksiter(0),
234  prev_max_k(0.),
235  prev_max_dzdl(0.),
236  prev_p_inv(0.),
237  seed_layer(0),
238  nthreads(num_threads),
239  vssp(NULL),
240  pins(NULL),
241  is_parallel(parallel),
242  is_thread(false),
243  ca_chi2_cut(2.0),
244  cosang_cut(0.985) {
245  vector<float> detector_material;
246 
247  for (unsigned int i = 0; i < radius.size(); ++i) {
248  detector_radii.push_back(radius[i]);
249  }
250  for (unsigned int i = 0; i < material.size(); ++i) {
251  detector_scatter.push_back(1.41421356237309515 * 0.0136 *
252  sqrt(3. * material[i]));
253  detector_material.push_back(3. * material[i]);
254  }
255 
256  detector_B_field = Bfield;
257 
258  integrated_scatter.assign(detector_scatter.size(), 0.);
259  float total_scatter_2 = 0.;
260  for (unsigned int l = 0; l < detector_scatter.size(); ++l) {
261  total_scatter_2 += detector_scatter[l] * detector_scatter[l];
262  integrated_scatter[l] = sqrt(total_scatter_2);
263  }
264 
265  kalman =
266  new CylinderKalman(detector_radii, detector_material, detector_B_field);
267 
268  vector<SimpleHit3D> one_layer;
269  layer_sorted.assign(n_layers, one_layer);
270  for (unsigned int i = 0; i < 4; ++i) {
271  layer_sorted_1[i].assign(n_layers, one_layer);
272  }
273  temp_comb.assign(n_layers, 0);
274 
275  if (is_parallel == true) {
276  Seamstress::init_vector(num_threads, vss);
277 
278  vssp = new vector<Seamstress*>();
279  for (unsigned int i = 0; i < vss.size(); i++) {
280  vssp->push_back(&(vss[i]));
281  }
282 
283  pins = new Pincushion<sPHENIXTracker>(this, vssp);
284 
285  vector<vector<unsigned int> > zoom_profile_new;
286  for (unsigned int i = 1; i < zoom_profile.size(); ++i) {
287  zoom_profile_new.push_back(zoom_profile[i]);
288  }
289 
290  for (unsigned int i = 0; i < nthreads; ++i) {
291  thread_trackers.push_back(new sPHENIXTracker(zoom_profile, minzoom, range,
292  material, radius, Bfield));
293  thread_trackers.back()->setThread();
294  thread_trackers.back()->setStartZoom(1);
295  thread_tracks.push_back(vector<SimpleTrack3D>());
296  thread_ranges.push_back(HelixRange());
297  thread_hits.push_back(vector<SimpleHit3D>());
298  split_output_hits.push_back(new vector<vector<SimpleHit3D> >());
299  split_ranges.push_back(new vector<HelixRange>());
300  split_input_hits.push_back(vector<SimpleHit3D>());
301  }
302  }
303 }
304 
306  if (kalman != NULL) delete kalman;
307  for (unsigned int i = 0; i < vss.size(); i++) {
308  vss[i].stop();
309  }
310  for (unsigned int i = 0; i < thread_trackers.size(); ++i) {
311  delete thread_trackers[i];
312  delete split_output_hits[i];
313  delete split_ranges[i];
314  }
315 
316  if (pins != NULL) delete pins;
317  if (vssp != NULL) delete vssp;
318 }
319 
320 float sPHENIXTracker::kappaToPt(float kappa) {
321  return detector_B_field / 333.6 / kappa;
322 }
323 
325  return detector_B_field / 333.6 / pt;
326 }
327 
328 // hel should be +- 1
329 // static void xyTangent(SimpleHit3D& hit1, SimpleHit3D& hit2, float kappa,
330 // float hel, float& ux_out, float& uy_out, float& ux_in,
331 // float& uy_in) {
332 // float x = hit2.get_x() - hit1.get_x();
333 // float y = hit2.get_y() - hit1.get_y();
334 // float D = sqrt(x * x + y * y);
335 // float ak = 0.5 * kappa * D;
336 // float D_inv = 1. / D;
337 // float hk = sqrt(1. - ak * ak);
338 
339 // float kcx = (ak * x + hel * hk * y) * D_inv;
340 // float kcy = (ak * y - hel * hk * x) * D_inv;
341 // float ktx = -(kappa * y - kcy);
342 // float kty = kappa * x - kcx;
343 // float norm = 1. / sqrt(ktx * ktx + kty * kty);
344 // ux_out = ktx * norm;
345 // uy_out = kty * norm;
346 
347 // ktx = kcy;
348 // kty = -kcx;
349 // norm = 1. / sqrt(ktx * ktx + kty * kty);
350 // ux_in = ktx * norm;
351 // uy_in = kty * norm;
352 // }
353 
354 // hel should be +- 1
355 // static float cosScatter(SimpleHit3D& hit1, SimpleHit3D& hit2, SimpleHit3D& hit3,
356 // float kappa, float hel) {
357 // float ux_in = 0.;
358 // float uy_in = 0.;
359 // float ux_out = 0.;
360 // float uy_out = 0.;
361 
362 // float temp1 = 0.;
363 // float temp2 = 0.;
364 
365 // xyTangent(hit1, hit2, kappa, hel, ux_in, uy_in, temp1, temp2);
366 // xyTangent(hit2, hit3, kappa, hel, temp1, temp2, ux_out, uy_out);
367 
368 // return ux_in * ux_out + uy_in * uy_out;
369 // }
370 
371 // static float dzdsSimple(SimpleHit3D& hit1, SimpleHit3D& hit2, float k) {
372 // float x = hit2.get_x() - hit1.get_x();
373 // float y = hit2.get_y() - hit1.get_y();
374 // float D = sqrt(x * x + y * y);
375 // float s = 0.;
376 // float temp1 = k * D * 0.5;
377 // temp1 *= temp1;
378 // float temp2 = D * 0.5;
379 // s += 2. * temp2;
380 // temp2 *= temp1;
381 // s += temp2 / 3.;
382 // temp2 *= temp1;
383 // s += (3. / 20.) * temp2;
384 // temp2 *= temp1;
385 // s += (5. / 56.) * temp2;
386 
387 // return (hit2.get_z() - hit1.get_z()) / s;
388 // }
389 
390 float sPHENIXTracker::dcaToVertexXY(SimpleTrack3D& track, float vx, float vy) {
391  float d_out = 0.;
392 
393  // find point at the dca to 0
394  float x0 = track.d * cos(track.phi);
395  float y0 = track.d * sin(track.phi);
396 
397  // change variables so x0,y0 -> 0,0
398  float phi2 =
399  atan2((1. + track.kappa * track.d) * sin(track.phi) - track.kappa * y0,
400  (1. + track.kappa * track.d) * cos(track.phi) - track.kappa * x0);
401 
402  // translate so that (0,0) -> (x0 - vx , y0 - vy)
403  float cosphi = cos(phi2);
404  float sinphi = sin(phi2);
405  float tx = cosphi + track.kappa * (x0 - vx);
406  float ty = sinphi + track.kappa * (y0 - vy);
407  float dk = sqrt(tx * tx + ty * ty) - 1.;
408  if (track.kappa == 0.) {
409  d_out = (x0 - vx) * cosphi + (y0 - vy) * sinphi;
410  } else {
411  d_out = dk / track.kappa;
412  }
413  return fabs(d_out);
414 }
415 
416 void sPHENIXTracker::finalize(vector<SimpleTrack3D>& input,
417  vector<SimpleTrack3D>& output) {
418 
419  if (is_thread == true) {
420  for (unsigned int i = 0; i < input.size(); ++i) {
421  output.push_back(input[i]);
422  }
423  return;
424  }
425 
426  unsigned int nt = input.size();
427  vector<bool> usetrack;
428  usetrack.assign(input.size(), true);
429  vector<float> next_best_chi2;
430  next_best_chi2.assign(input.size(), 99999.);
431 
432  if (reject_ghosts == true) {
433  tripletRejection(input, output, usetrack, next_best_chi2);
434  }
435 
436  vector<HelixKalmanState> states_new;
437 
438  for (unsigned int i = 0; i < nt; ++i) {
439  if (usetrack[i] == true) {
440  if (!(track_states[i].chi2 == track_states[i].chi2)) {
441  continue;
442  }
443 
444  output.push_back(input[i]);
445  output.back().index = (output.size() - 1);
446  states_new.push_back(track_states[i]);
447  isolation_variable.push_back(next_best_chi2[i]);
448  }
449  }
450  track_states = states_new;
451  if (smooth_back == true) {
452  for (unsigned int i = 0; i < output.size(); ++i) {
453 
454  HelixKalmanState state = track_states[i];
455 
456  track_states[i].C *= 3.;
457  track_states[i].chi2 = 0.;
458  track_states[i].x_int = 0.;
459  track_states[i].y_int = 0.;
460  track_states[i].z_int = 0.;
461  // track_states[i].position = output[i].hits.size();
462  track_states[i].position = 0;
463  for (int h = (output[i].hits.size() - 1); h >= 0; --h) {
464  SimpleHit3D hit = output[i].hits[h];
465  float err_scale = 1.;
466  int layer = hit.get_layer();
467  if ((layer >= 0) && (layer < (int)(hit_error_scale.size()))) {
468  err_scale = hit_error_scale[layer];
469  }
470  err_scale *=
471  3.0; // fudge factor, like needed due to non-gaussian errors
472 
473  // \todo location of a rescale fudge factor
474 
475  kalman->addHit(hit, track_states[i]);
476  track_states[i].position = h;
477  }
478 
479  // track_states[i].C = state.C;
480  track_states[i].chi2 = state.chi2;
481  track_states[i].C *= 2. / 3.;
482 
483  output[i].phi = track_states[i].phi;
484  output[i].d = track_states[i].d;
485  output[i].kappa = track_states[i].kappa;
486  output[i].z0 = track_states[i].z0;
487  output[i].dzdl = track_states[i].dzdl;
488  }
489  }
490 
491  if (verbosity > 0) {
492  cout << "# fits = " << nfits << endl;
493  cout << "findTracks called " << findtracksiter << " times" << endl;
494  cout << "CAtime = " << CAtime << endl;
495  cout << "KALime = " << KALtime << endl;
496  }
497 }
498 
499 void sPHENIXTracker::findTracks(vector<SimpleHit3D>& hits,
500  vector<SimpleTrack3D>& tracks,
501  const HelixRange& range) {
502  findtracksiter += 1;
503  findTracksBySegments(hits, tracks, range);
504 }
505 
506 bool sPHENIXTracker::breakRecursion(const vector<SimpleHit3D>& hits,
507  const HelixRange& /*range*/) {
508  if (seeding == true) {
509  return false;
510  }
511  unsigned int layer_mask[4] = {0, 0, 0, 0};
512  for (unsigned int i = 0; i < hits.size(); ++i) {
513  if (hits[i].get_layer() < 32) {
514  layer_mask[0] = layer_mask[0] | (1 << hits[i].get_layer());
515  } else if (hits[i].get_layer() < 64) {
516  layer_mask[1] = layer_mask[1] | (1 << (hits[i].get_layer() - 32));
517  } else if (hits[i].get_layer() < 96) {
518  layer_mask[2] = layer_mask[2] | (1 << (hits[i].get_layer() - 64));
519  } else if (hits[i].get_layer() < 128) {
520  layer_mask[3] = layer_mask[3] | (1 << (hits[i].get_layer() - 96));
521  }
522  }
523  unsigned int nlayers =
524  __builtin_popcount(layer_mask[0]) + __builtin_popcount(layer_mask[1]) +
525  __builtin_popcount(layer_mask[2]) + __builtin_popcount(layer_mask[3]);
526 
527  return (nlayers < required_layers);
528 }
529 
530 float sPHENIXTracker::phiError(SimpleHit3D& hit, float /*min_k*/, float max_k,
531  float min_d, float max_d, float /*min_z0*/,
532  float /*max_z0*/, float /*min_dzdl*/, float max_dzdl,
533  bool pairvoting) {
534  float Bfield_inv = 1. / detector_B_field;
535  float p_inv = 0.;
536 
537  if ((prev_max_k == max_k) && (prev_max_dzdl == max_dzdl)) {
538  p_inv = prev_p_inv;
539  } else {
540  prev_max_k = max_k;
541  prev_max_dzdl = max_dzdl;
542  prev_p_inv = 3.33333333333333314e+02 * max_k * Bfield_inv *
543  sqrt(1. - max_dzdl * max_dzdl);
544  p_inv = prev_p_inv;
545  }
546  float total_scatter_2 = 0.;
547  for (int i = seed_layer + 1; i <= (hit.get_layer()); ++i) {
548  float this_scatter = detector_scatter[i - 1] *
549  (detector_radii[i] - detector_radii[i - 1]) /
550  detector_radii[i];
551  total_scatter_2 += this_scatter * this_scatter;
552  }
553  float angle = p_inv * sqrt(total_scatter_2) * 1.0;
554  float dsize = 0.5 * (max_d - min_d);
555  float angle_from_d = dsize / detector_radii[hit.get_layer()];
556  float returnval = 0.;
557  if (pairvoting == false) {
558  if (angle_from_d > angle) {
559  returnval = 0.;
560  } else {
561  returnval = (angle - angle_from_d);
562  }
563  } else {
564  returnval = angle;
565  }
566 
567  return returnval;
568 }
569 
570 float sPHENIXTracker::dzdlError(SimpleHit3D& hit, float /*min_k*/, float max_k,
571  float /*min_d*/, float /*max_d*/, float min_z0,
572  float max_z0, float /*min_dzdl*/, float max_dzdl,
573  bool pairvoting) {
574  float Bfield_inv = 1. / detector_B_field;
575  float p_inv = 0.;
576 
577  if ((prev_max_k == max_k) && (prev_max_dzdl == max_dzdl)) {
578  p_inv = prev_p_inv;
579  } else {
580  prev_max_k = max_k;
581  prev_max_dzdl = max_dzdl;
582  prev_p_inv = 3.33333333333333314e+02 * max_k * Bfield_inv *
583  sqrt(1. - max_dzdl * max_dzdl);
584  p_inv = prev_p_inv;
585  }
586  float total_scatter_2 = 0.;
587  for (int i = seed_layer + 1; i <= (hit.get_layer()); ++i) {
588  float this_scatter = detector_scatter[i - 1] *
589  (detector_radii[i] - detector_radii[i - 1]) /
590  detector_radii[i];
591  total_scatter_2 += this_scatter * this_scatter;
592  }
593  float angle = p_inv * sqrt(total_scatter_2) * 1.0;
594  float z0size = 0.5 * (max_z0 - min_z0);
595  float angle_from_z0 = z0size / detector_radii[hit.get_layer()];
596  float returnval = 0.;
597  if (pairvoting == false) {
598  if (angle_from_z0 > angle) {
599  returnval = 0.;
600  } else {
601  returnval = (angle - angle_from_z0);
602  }
603  } else {
604  returnval = angle;
605  }
606 
607  return returnval;
608 }
609 
611  HelixKalmanState* state = &(seed_states[seed.index]);
612 
613  float dphi = 2. * sqrt(state->C(0, 0));
614  float dd = 2. * sqrt(state->C(1, 1));
615  float dk = 2. * state->C(2, 2);
616  float dz0 = 2. * sqrt(state->C(3, 3));
617  float ddzdl = 2. * sqrt(state->C(4, 4));
618 
619  range.min_phi = seed.phi - dphi;
620  range.max_phi = seed.phi + dphi;
621  if (range.min_phi < 0.) {
622  range.min_phi = 0.;
623  }
624  if (range.max_phi > 2. * M_PI) {
625  range.max_phi = 2. * M_PI;
626  }
627  range.min_d = seed.d - dd;
628  range.max_d = seed.d + dd;
629  range.min_k = seed.kappa - dk;
630  range.max_k = seed.kappa + dk;
631  if (range.min_k < 0.) {
632  range.min_k = 0.;
633  }
634 
635  range.min_k = range.min_k * range.min_k;
636  range.max_k = range.max_k * range.max_k;
637 
638  range.min_dzdl = seed.dzdl - ddzdl;
639  range.max_dzdl = seed.dzdl + ddzdl;
640  range.min_z0 = seed.z0 - dz0;
641  range.max_z0 = seed.z0 + dz0;
642 }
643 
644 
646  vector<float> chi2_hit;
647  return sPHENIXTracker::fitTrack(track, chi2_hit);
648 }
649 
650 float sPHENIXTracker::fitTrack(SimpleTrack3D& track, vector<float>& chi2_hit) {
651  chi2_hit.clear();
652  vector<float> xyres;
653  vector<float> xyres_inv;
654  vector<float> zres;
655  vector<float> zres_inv;
656  for (unsigned int i = 0; i < track.hits.size(); i++) {
657 
659 
660  xyres.push_back(sqrt((0.5*sqrt(12.)*sqrt(track.hits[i].get_size(0,0))) * (0.5*sqrt(12.)*sqrt(track.hits[i].get_size(0,0))) +
661  (0.5*sqrt(12.)*sqrt(track.hits[i].get_size(1,1))) * (0.5*sqrt(12.)*sqrt(track.hits[i].get_size(1,1)))));
662  xyres_inv.push_back(1. / xyres.back());
663  zres.push_back((0.5*sqrt(12.)*sqrt(track.hits[i].get_size(2,2))));
664  zres_inv.push_back(1. / zres.back());
665  }
666 
667  chi2_hit.resize(track.hits.size(), 0.);
668 
669  MatrixXf y = MatrixXf::Zero(track.hits.size(), 1);
670  for (unsigned int i = 0; i < track.hits.size(); i++) {
671  y(i, 0) = (pow(track.hits[i].get_x(), 2) + pow(track.hits[i].get_y(), 2));
672  y(i, 0) *= xyres_inv[i];
673  }
674 
675  MatrixXf X = MatrixXf::Zero(track.hits.size(), 3);
676  for (unsigned int i = 0; i < track.hits.size(); i++) {
677  X(i, 0) = track.hits[i].get_x();
678  X(i, 1) = track.hits[i].get_y();
679  X(i, 2) = -1.;
680  X(i, 0) *= xyres_inv[i];
681  X(i, 1) *= xyres_inv[i];
682  X(i, 2) *= xyres_inv[i];
683  }
684 
685  MatrixXf Xt = X.transpose();
686 
687  MatrixXf prod = Xt * X;
688 
689  MatrixXf Xty = Xt * y;
690  MatrixXf beta = prod.ldlt().solve(Xty);
691 
692  float cx = beta(0, 0) * 0.5;
693  float cy = beta(1, 0) * 0.5;
694  float r = sqrt(cx * cx + cy * cy - beta(2, 0));
695 
696  float phi = atan2(cy, cx);
697  float d = sqrt(cx * cx + cy * cy) - r;
698  float k = 1. / r;
699 
700  MatrixXf diff = y - (X * beta);
701  MatrixXf chi2 = (diff.transpose()) * diff;
702 
703  float dx = d * cos(phi);
704  float dy = d * sin(phi);
705 
706  MatrixXf y2 = MatrixXf::Zero(track.hits.size(), 1);
707  for (unsigned int i = 0; i < track.hits.size(); i++) {
708  y2(i, 0) = track.hits[i].get_z();
709  y2(i, 0) *= zres_inv[i];
710  }
711 
712  MatrixXf X2 = MatrixXf::Zero(track.hits.size(), 2);
713  for (unsigned int i = 0; i < track.hits.size(); i++) {
714  float D = sqrt(pow(dx - track.hits[i].get_x(), 2) + pow(dy - track.hits[i].get_y(), 2));
715  float s = 0.0;
716 
717  if (0.5 * k * D > 0.1) {
718  float v = 0.5 * k * D;
719  if (v >= 0.999999) {
720  v = 0.999999;
721  }
722  s = 2. * asin(v) / k;
723  } else {
724  float temp1 = k * D * 0.5;
725  temp1 *= temp1;
726  float temp2 = D * 0.5;
727  s += 2. * temp2;
728  temp2 *= temp1;
729  s += temp2 / 3.;
730  temp2 *= temp1;
731  s += (3. / 20.) * temp2;
732  temp2 *= temp1;
733  s += (5. / 56.) * temp2;
734  }
735 
736  X2(i, 0) = s;
737  X2(i, 1) = 1.0;
738 
739  X2(i, 0) *= zres_inv[i];
740  X2(i, 1) *= zres_inv[i];
741  }
742 
743  MatrixXf Xt2 = X2.transpose();
744  MatrixXf prod2 = Xt2 * X2;
745 
746  MatrixXf Xty2 = Xt2 * y2;
747  MatrixXf beta2 = prod2.ldlt().solve(Xty2);
748 
749  MatrixXf diff2 = y2 - (X2 * beta2);
750  MatrixXf chi2_z = (diff2.transpose()) * diff2;
751 
752  float z0 = beta2(1, 0);
753  float dzdl = beta2(0, 0) / sqrt(1. + beta2(0, 0) * beta2(0, 0));
754 
755  track.phi = phi;
756  track.d = d;
757  track.kappa = k;
758  track.dzdl = dzdl;
759  track.z0 = z0;
760 
761  if (track.kappa != 0.) {
762  r = 1. / track.kappa;
763  } else {
764  r = 1.0e10;
765  }
766 
767  cx = (track.d + r) * cos(track.phi);
768  cy = (track.d + r) * sin(track.phi);
769 
770  float chi2_tot = 0.;
771  for (unsigned int h = 0; h < track.hits.size(); h++) {
772  float dx1 = track.hits[h].get_x() - cx;
773  float dy1 = track.hits[h].get_y() - cy;
774 
775  float dx2 = track.hits[h].get_x() + cx;
776  float dy2 = track.hits[h].get_y() + cy;
777 
778  float xydiff1 = sqrt(dx1 * dx1 + dy1 * dy1) - r;
779  float xydiff2 = sqrt(dx2 * dx2 + dy2 * dy2) - r;
780  float xydiff = xydiff2;
781  if (fabs(xydiff1) < fabs(xydiff2)) {
782  xydiff = xydiff1;
783  }
784 
785  float ls_xy = xyres[h];
786 
787  chi2_hit[h] = 0.;
788  chi2_hit[h] += xydiff * xydiff / (ls_xy * ls_xy);
789  chi2_hit[h] += diff2(h, 0) * diff2(h, 0);
790 
791  chi2_tot += chi2_hit[h];
792  }
793 
794  unsigned int deg_of_freedom = 2 * track.hits.size() - 5;
795 
796  return (chi2_tot) / ((float)(deg_of_freedom));
797 }
798 
800  float* x1_a, float* y1_a, float* z1_a, float* x2_a, float* y2_a,
801  float* z2_a, float* x3_a, float* y3_a, float* z3_a, float* dx1_a,
802  float* dy1_a, float* dz1_a, float* dx2_a, float* dy2_a, float* dz2_a,
803  float* dx3_a, float* dy3_a, float* dz3_a, float* kappa_a, float* dkappa_a,
804  float* ux_mid_a, float* uy_mid_a, float* ux_end_a, float* uy_end_a,
805  float* dzdl_1_a, float* dzdl_2_a, float* ddzdl_1_a, float* ddzdl_2_a) {
806  static const __m128 two = {2., 2., 2., 2.};
807 
808  __m128 x1 = _mm_load_ps(x1_a);
809  __m128 x2 = _mm_load_ps(x2_a);
810  __m128 x3 = _mm_load_ps(x3_a);
811  __m128 y1 = _mm_load_ps(y1_a);
812  __m128 y2 = _mm_load_ps(y2_a);
813  __m128 y3 = _mm_load_ps(y3_a);
814  __m128 z1 = _mm_load_ps(z1_a);
815  __m128 z2 = _mm_load_ps(z2_a);
816  __m128 z3 = _mm_load_ps(z3_a);
817 
818  __m128 dx1 = _mm_load_ps(dx1_a);
819  __m128 dx2 = _mm_load_ps(dx2_a);
820  __m128 dx3 = _mm_load_ps(dx3_a);
821  __m128 dy1 = _mm_load_ps(dy1_a);
822  __m128 dy2 = _mm_load_ps(dy2_a);
823  __m128 dy3 = _mm_load_ps(dy3_a);
824  __m128 dz1 = _mm_load_ps(dz1_a);
825  __m128 dz2 = _mm_load_ps(dz2_a);
826  __m128 dz3 = _mm_load_ps(dz3_a);
827 
828  __m128 D12 = _mm_sub_ps(x2, x1);
829  D12 = _mm_mul_ps(D12, D12);
830  __m128 tmp1 = _mm_sub_ps(y2, y1);
831  tmp1 = _mm_mul_ps(tmp1, tmp1);
832  D12 = _mm_add_ps(D12, tmp1);
833  D12 = _vec_sqrt_ps(D12);
834 
835  __m128 D23 = _mm_sub_ps(x3, x2);
836  D23 = _mm_mul_ps(D23, D23);
837  tmp1 = _mm_sub_ps(y3, y2);
838  tmp1 = _mm_mul_ps(tmp1, tmp1);
839  D23 = _mm_add_ps(D23, tmp1);
840  D23 = _vec_sqrt_ps(D23);
841 
842  __m128 D31 = _mm_sub_ps(x1, x3);
843  D31 = _mm_mul_ps(D31, D31);
844  tmp1 = _mm_sub_ps(y1, y3);
845  tmp1 = _mm_mul_ps(tmp1, tmp1);
846  D31 = _mm_add_ps(D31, tmp1);
847  D31 = _vec_sqrt_ps(D31);
848 
849  __m128 k = _mm_mul_ps(D12, D23);
850  k = _mm_mul_ps(k, D31);
851  k = _vec_rec_ps(k);
852  tmp1 = (D12 + D23 + D31) * (D23 + D31 - D12) * (D12 + D31 - D23) *
853  (D12 + D23 - D31);
854  tmp1 = _vec_sqrt_ps(tmp1);
855  k *= tmp1;
856 
857  __m128 tmp2 = _mm_cmpgt_ps(tmp1, zero);
858  tmp1 = _mm_and_ps(tmp2, k);
859  tmp2 = _mm_andnot_ps(tmp2, zero);
860  k = _mm_xor_ps(tmp1, tmp2);
861 
862  _mm_store_ps(kappa_a, k);
863  __m128 k_inv = _vec_rec_ps(k);
864 
865  __m128 D12_inv = _vec_rec_ps(D12);
866  __m128 D23_inv = _vec_rec_ps(D23);
867  __m128 D31_inv = _vec_rec_ps(D31);
868 
869  __m128 dr1 = dx1 * dx1 + dy1 * dy1;
870  dr1 = _vec_sqrt_ps(dr1);
871  __m128 dr2 = dx2 * dx2 + dy2 * dy2;
872  dr2 = _vec_sqrt_ps(dr2);
873  __m128 dr3 = dx3 * dx3 + dy3 * dy3;
874  dr3 = _vec_sqrt_ps(dr3);
875 
876  __m128 dk1 = (dr1 + dr2) * D12_inv * D12_inv;
877  __m128 dk2 = (dr2 + dr3) * D23_inv * D23_inv;
878  __m128 dk = dk1 + dk2;
879  _mm_store_ps(dkappa_a, dk);
880 
881  __m128 ux12 = (x2 - x1) * D12_inv;
882  __m128 uy12 = (y2 - y1) * D12_inv;
883  __m128 ux23 = (x3 - x2) * D23_inv;
884  __m128 uy23 = (y3 - y2) * D23_inv;
885  __m128 ux13 = (x3 - x1) * D31_inv;
886  __m128 uy13 = (y3 - y1) * D31_inv;
887 
888  __m128 cosalpha = ux12 * ux13 + uy12 * uy13;
889  __m128 sinalpha = ux13 * uy12 - ux12 * uy13;
890 
891  __m128 ux_mid = ux23 * cosalpha - uy23 * sinalpha;
892  __m128 uy_mid = ux23 * sinalpha + uy23 * cosalpha;
893  _mm_store_ps(ux_mid_a, ux_mid);
894  _mm_store_ps(uy_mid_a, uy_mid);
895 
896  __m128 ux_end = ux23 * cosalpha + uy23 * sinalpha;
897  __m128 uy_end = uy23 * cosalpha - ux23 * sinalpha;
898 
899  _mm_store_ps(ux_end_a, ux_end);
900  _mm_store_ps(uy_end_a, uy_end);
901 
902  // asin(x) = 2*atan( x/( 1 + sqrt( 1 - x*x ) ) )
903  __m128 v = one - sinalpha * sinalpha;
904  v = _vec_sqrt_ps(v);
905  v += one;
906  v = _vec_rec_ps(v);
907  v *= sinalpha;
908  __m128 s2 = _vec_atan_ps(v);
909  s2 *= two;
910  s2 *= k_inv;
911  tmp1 = _mm_cmpgt_ps(k, zero);
912  tmp2 = _mm_and_ps(tmp1, s2);
913  tmp1 = _mm_andnot_ps(tmp1, D23);
914  s2 = _mm_xor_ps(tmp1, tmp2);
915 
916  // dz/dl = (dz/ds)/sqrt(1 + (dz/ds)^2)
917  // = dz/sqrt(s^2 + dz^2)
918  __m128 del_z_2 = z3 - z2;
919  __m128 dzdl_2 = s2 * s2 + del_z_2 * del_z_2;
920  dzdl_2 = _vec_rsqrt_ps(dzdl_2);
921  dzdl_2 *= del_z_2;
922  __m128 ddzdl_2 = (dz2 + dz3) * D23_inv;
923  _mm_store_ps(dzdl_2_a, dzdl_2);
924  _mm_store_ps(ddzdl_2_a, ddzdl_2);
925 
926  sinalpha = ux13 * uy23 - ux23 * uy13;
927  v = one - sinalpha * sinalpha;
928  v = _vec_sqrt_ps(v);
929  v += one;
930  v = _vec_rec_ps(v);
931  v *= sinalpha;
932  __m128 s1 = _vec_atan_ps(v);
933  s1 *= two;
934  s1 *= k_inv;
935  tmp1 = _mm_cmpgt_ps(k, zero);
936  tmp2 = _mm_and_ps(tmp1, s1);
937  tmp1 = _mm_andnot_ps(tmp1, D12);
938  s1 = _mm_xor_ps(tmp1, tmp2);
939 
940  __m128 del_z_1 = z2 - z1;
941  __m128 dzdl_1 = s1 * s1 + del_z_1 * del_z_1;
942  dzdl_1 = _vec_rsqrt_ps(dzdl_1);
943  dzdl_1 *= del_z_1;
944  __m128 ddzdl_1 = (dz1 + dz2) * D12_inv;
945  _mm_store_ps(dzdl_1_a, dzdl_1);
946  _mm_store_ps(ddzdl_1_a, ddzdl_1);
947 }
948 
950  float* x1_a, float* y1_a, float* z1_a, float* x2_a, float* y2_a,
951  float* z2_a, float* x3_a, float* y3_a, float* z3_a, float* dx1_a,
952  float* dy1_a, float* dz1_a, float* dx2_a, float* dy2_a, float* dz2_a,
953  float* dx3_a, float* dy3_a, float* dz3_a, float* kappa_a, float* dkappa_a,
954  float* ux_mid_a, float* uy_mid_a, float* ux_end_a, float* uy_end_a,
955  float* dzdl_1_a, float* dzdl_2_a, float* ddzdl_1_a, float* ddzdl_2_a,
956  float sinang_cut, float cosang_diff_inv, float* cur_kappa_a,
957  float* cur_dkappa_a, float* cur_ux_a, float* cur_uy_a, float* cur_chi2_a,
958  float* chi2_a) {
959  static const __m128 two = {2., 2., 2., 2.};
960 
961  __m128 x1 = _mm_load_ps(x1_a);
962  __m128 x2 = _mm_load_ps(x2_a);
963  __m128 x3 = _mm_load_ps(x3_a);
964  __m128 y1 = _mm_load_ps(y1_a);
965  __m128 y2 = _mm_load_ps(y2_a);
966  __m128 y3 = _mm_load_ps(y3_a);
967  __m128 z1 = _mm_load_ps(z1_a);
968  __m128 z2 = _mm_load_ps(z2_a);
969  __m128 z3 = _mm_load_ps(z3_a);
970 
971  __m128 dx1 = _mm_load_ps(dx1_a);
972  __m128 dx2 = _mm_load_ps(dx2_a);
973  __m128 dx3 = _mm_load_ps(dx3_a);
974  __m128 dy1 = _mm_load_ps(dy1_a);
975  __m128 dy2 = _mm_load_ps(dy2_a);
976  __m128 dy3 = _mm_load_ps(dy3_a);
977  __m128 dz1 = _mm_load_ps(dz1_a);
978  __m128 dz2 = _mm_load_ps(dz2_a);
979  __m128 dz3 = _mm_load_ps(dz3_a);
980 
981  __m128 D12 = _mm_sub_ps(x2, x1);
982  D12 = _mm_mul_ps(D12, D12);
983  __m128 tmp1 = _mm_sub_ps(y2, y1);
984  tmp1 = _mm_mul_ps(tmp1, tmp1);
985  D12 = _mm_add_ps(D12, tmp1);
986  D12 = _vec_sqrt_ps(D12);
987 
988  __m128 D23 = _mm_sub_ps(x3, x2);
989  D23 = _mm_mul_ps(D23, D23);
990  tmp1 = _mm_sub_ps(y3, y2);
991  tmp1 = _mm_mul_ps(tmp1, tmp1);
992  D23 = _mm_add_ps(D23, tmp1);
993  D23 = _vec_sqrt_ps(D23);
994 
995  __m128 D31 = _mm_sub_ps(x1, x3);
996  D31 = _mm_mul_ps(D31, D31);
997  tmp1 = _mm_sub_ps(y1, y3);
998  tmp1 = _mm_mul_ps(tmp1, tmp1);
999  D31 = _mm_add_ps(D31, tmp1);
1000  D31 = _vec_sqrt_ps(D31);
1001 
1002  __m128 k = _mm_mul_ps(D12, D23);
1003  k = _mm_mul_ps(k, D31);
1004  k = _vec_rec_ps(k);
1005  tmp1 = (D12 + D23 + D31) * (D23 + D31 - D12) * (D12 + D31 - D23) *
1006  (D12 + D23 - D31);
1007  tmp1 = _vec_sqrt_ps(tmp1);
1008  k *= tmp1;
1009 
1010  __m128 tmp2 = _mm_cmpgt_ps(tmp1, zero);
1011  tmp1 = _mm_and_ps(tmp2, k);
1012  tmp2 = _mm_andnot_ps(tmp2, zero);
1013  k = _mm_xor_ps(tmp1, tmp2);
1014 
1015  _mm_store_ps(kappa_a, k);
1016  __m128 k_inv = _vec_rec_ps(k);
1017 
1018  __m128 D12_inv = _vec_rec_ps(D12);
1019  __m128 D23_inv = _vec_rec_ps(D23);
1020  __m128 D31_inv = _vec_rec_ps(D31);
1021 
1022  __m128 dr1 = dx1 * dx1 + dy1 * dy1;
1023  dr1 = _vec_sqrt_ps(dr1);
1024  __m128 dr2 = dx2 * dx2 + dy2 * dy2;
1025  dr2 = _vec_sqrt_ps(dr2);
1026  __m128 dr3 = dx3 * dx3 + dy3 * dy3;
1027  dr3 = _vec_sqrt_ps(dr3);
1028 
1029  __m128 dk1 = (dr1 + dr2) * D12_inv * D12_inv;
1030  __m128 dk2 = (dr2 + dr3) * D23_inv * D23_inv;
1031  __m128 dk = dk1 + dk2;
1032  _mm_store_ps(dkappa_a, dk);
1033 
1034  __m128 ux12 = (x2 - x1) * D12_inv;
1035  __m128 uy12 = (y2 - y1) * D12_inv;
1036  __m128 ux23 = (x3 - x2) * D23_inv;
1037  __m128 uy23 = (y3 - y2) * D23_inv;
1038  __m128 ux13 = (x3 - x1) * D31_inv;
1039  __m128 uy13 = (y3 - y1) * D31_inv;
1040 
1041  __m128 cosalpha = ux12 * ux13 + uy12 * uy13;
1042  __m128 sinalpha = ux13 * uy12 - ux12 * uy13;
1043 
1044  __m128 ux_mid = ux23 * cosalpha - uy23 * sinalpha;
1045  __m128 uy_mid = ux23 * sinalpha + uy23 * cosalpha;
1046  _mm_store_ps(ux_mid_a, ux_mid);
1047  _mm_store_ps(uy_mid_a, uy_mid);
1048 
1049  __m128 ux_end = ux23 * cosalpha + uy23 * sinalpha;
1050  __m128 uy_end = uy23 * cosalpha - ux23 * sinalpha;
1051 
1052  _mm_store_ps(ux_end_a, ux_end);
1053  _mm_store_ps(uy_end_a, uy_end);
1054 
1055  // asin(x) = 2*atan( x/( 1 + sqrt( 1 - x*x ) ) )
1056  __m128 v = one - sinalpha * sinalpha;
1057  v = _vec_sqrt_ps(v);
1058  v += one;
1059  v = _vec_rec_ps(v);
1060  v *= sinalpha;
1061  __m128 s2 = _vec_atan_ps(v);
1062  s2 *= two;
1063  s2 *= k_inv;
1064  tmp1 = _mm_cmpgt_ps(k, zero);
1065  tmp2 = _mm_and_ps(tmp1, s2);
1066  tmp1 = _mm_andnot_ps(tmp1, D23);
1067  s2 = _mm_xor_ps(tmp1, tmp2);
1068 
1069  // dz/dl = (dz/ds)/sqrt(1 + (dz/ds)^2)
1070  // = dz/sqrt(s^2 + dz^2)
1071  __m128 del_z_2 = z3 - z2;
1072  __m128 dzdl_2 = s2 * s2 + del_z_2 * del_z_2;
1073  dzdl_2 = _vec_rsqrt_ps(dzdl_2);
1074  dzdl_2 *= del_z_2;
1075  __m128 ddzdl_2 = (dz2 + dz3) * D23_inv;
1076  _mm_store_ps(dzdl_2_a, dzdl_2);
1077  _mm_store_ps(ddzdl_2_a, ddzdl_2);
1078 
1079  sinalpha = ux13 * uy23 - ux23 * uy13;
1080  v = one - sinalpha * sinalpha;
1081  v = _vec_sqrt_ps(v);
1082  v += one;
1083  v = _vec_rec_ps(v);
1084  v *= sinalpha;
1085  __m128 s1 = _vec_atan_ps(v);
1086  s1 *= two;
1087  s1 *= k_inv;
1088  tmp1 = _mm_cmpgt_ps(k, zero);
1089  tmp2 = _mm_and_ps(tmp1, s1);
1090  tmp1 = _mm_andnot_ps(tmp1, D12);
1091  s1 = _mm_xor_ps(tmp1, tmp2);
1092 
1093  __m128 del_z_1 = z2 - z1;
1094  __m128 dzdl_1 = s1 * s1 + del_z_1 * del_z_1;
1095  dzdl_1 = _vec_rsqrt_ps(dzdl_1);
1096  dzdl_1 *= del_z_1;
1097  __m128 ddzdl_1 = (dz1 + dz2) * D12_inv;
1098  _mm_store_ps(dzdl_1_a, dzdl_1);
1099  _mm_store_ps(ddzdl_1_a, ddzdl_1);
1100 
1101  __m128 c_dk = _mm_load_ps(cur_dkappa_a);
1102  __m128 c_k = _mm_load_ps(cur_kappa_a);
1103  __m128 c_ux = _mm_load_ps(cur_ux_a);
1104  __m128 c_uy = _mm_load_ps(cur_uy_a);
1105  __m128 c_chi2 = _mm_load_ps(cur_chi2_a);
1106  __m128 sinang = _mm_load1_ps(&sinang_cut);
1107  __m128 cosdiff = _mm_load1_ps(&cosang_diff_inv);
1108 
1109  __m128 kdiff = c_k - k;
1110  __m128 n_dk = c_dk + dk + sinang * k;
1111  __m128 chi2_k = kdiff * kdiff / (n_dk * n_dk);
1112  __m128 cos_scatter = c_ux * ux_mid + c_uy * uy_mid;
1113  __m128 chi2_ang =
1114  (one - cos_scatter) * (one - cos_scatter) * cosdiff * cosdiff;
1115  tmp1 = dzdl_1 * sinang;
1116  _vec_fabs_ps(tmp1);
1117  __m128 chi2_dzdl = (dzdl_1 - dzdl_2) / (ddzdl_1 + ddzdl_2 + tmp1);
1118  chi2_dzdl *= chi2_dzdl;
1119  chi2_dzdl *= one_o_2;
1120 
1121  __m128 n_chi2 = c_chi2 + chi2_ang + chi2_k + chi2_dzdl;
1122  _mm_store_ps(chi2_a, n_chi2);
1123 }
1124 
1125 struct TempComb {
1129 
1130  inline bool operator<(const TempComb& other) const {
1131  if (track.hits.size() > other.track.hits.size()) {
1132  return true;
1133  } else if (track.hits.size() == other.track.hits.size()) {
1134  return state.chi2 < other.state.chi2;
1135  } else {
1136  return false;
1137  }
1138  }
1139 };
1140 
1141 void sPHENIXTracker::initDummyHits(vector<SimpleHit3D>& dummies,
1142  const HelixRange& range,
1143  HelixKalmanState& init_state) {
1144  SimpleTrack3D dummy_track;
1145  dummy_track.hits.push_back(SimpleHit3D());
1146  dummy_track.kappa = 0.5 * (range.min_k + range.max_k);
1147  dummy_track.phi = 0.5 * (range.min_phi + range.max_phi);
1148  dummy_track.d = 0.5 * (range.min_d + range.max_d);
1149  dummy_track.dzdl = 0.5 * (range.min_dzdl + range.max_dzdl);
1150  dummy_track.z0 = 0.5 * (range.min_z0 + range.max_z0);
1151 
1152  init_state.kappa = dummy_track.kappa;
1153  init_state.nu = sqrt(dummy_track.kappa);
1154  init_state.phi = dummy_track.phi;
1155  init_state.d = dummy_track.d;
1156  init_state.dzdl = dummy_track.dzdl;
1157  init_state.z0 = dummy_track.z0;
1158 
1159  init_state.C = Matrix<float, 5, 5>::Zero(5, 5);
1160  init_state.C(0, 0) = pow(range.max_phi - range.min_phi, 2.);
1161  init_state.C(1, 1) = pow(range.max_d - range.min_d, 2.);
1162  init_state.C(2, 2) = pow(10. * sqrt(range.max_k - range.min_k), 2.);
1163  init_state.C(3, 3) = pow(range.max_z0 - range.min_z0, 2.);
1164  init_state.C(4, 4) = pow(range.max_dzdl - range.min_dzdl, 2.);
1165  init_state.chi2 = 0.;
1166  init_state.position = 0;
1167  init_state.x_int = 0.;
1168  init_state.y_int = 0.;
1169  init_state.z_int = 0.;
1170 
1171  for (unsigned int i = 0; i < n_layers; ++i) {
1172  float x, y, z;
1173  projectToLayer(dummy_track, i, x, y, z);
1174  dummies[i].set_x(x);
1175  dummies[i].set_y(x);
1176  dummies[i].set_z(x);
1177  dummies[i].set_layer(i);
1178  }
1179 }
1180 
1181 
1182 
1183 // static void triplet_rejection(vector<SimpleTrack3D>& input,
1184 // vector<float>& chi2s, vector<bool>& usetrack) {
1185 // vector<hit_triplet> trips;
1186 // for (unsigned int i = 0; i < input.size(); ++i) {
1187 // for (unsigned int h1 = 0; h1 < input[i].hits.size(); ++h1) {
1188 // for (unsigned int h2 = (h1 + 1); h2 < input[i].hits.size(); ++h2) {
1189 // for (unsigned int h3 = (h2 + 1); h3 < input[i].hits.size(); ++h3) {
1190 // trips.push_back(hit_triplet(input[i].hits[h1].get_id(),
1191 // input[i].hits[h2].get_id(),
1192 // input[i].hits[h3].get_id(), i, chi2s[i]));
1193 // }
1194 // }
1195 // }
1196 // }
1197 // if (trips.size() == 0) {
1198 // return;
1199 // }
1200 // sort(trips.begin(), trips.end());
1201 // unsigned int pos = 0;
1202 // unsigned int cur_h1 = trips[pos].hit1;
1203 // unsigned int cur_h2 = trips[pos].hit2;
1204 // while (pos < trips.size()) {
1205 // unsigned int next_pos = pos + 1;
1206 // if (next_pos >= trips.size()) {
1207 // break;
1208 // }
1209 // while (trips[pos] == trips[next_pos]) {
1210 // next_pos += 1;
1211 // if (next_pos >= trips.size()) {
1212 // break;
1213 // }
1214 // }
1215 // if ((next_pos - pos) > 1) {
1216 // float best_chi2 = trips[pos].chi2;
1217 // float next_chi2 = trips[pos + 1].chi2;
1218 // unsigned int best_pos = pos;
1219 // for (unsigned int i = (pos + 1); i < next_pos; ++i) {
1220 // if (input[trips[i].track].hits.size() <
1221 // input[trips[best_pos].track].hits.size()) {
1222 // continue;
1223 // } else if ((input[trips[i].track].hits.size() >
1224 // input[trips[best_pos].track].hits.size()) ||
1225 // (input[trips[i].track].hits.back().get_layer() >
1226 // input[trips[best_pos].track].hits.back().get_layer())) {
1227 // next_chi2 = best_chi2;
1228 // best_chi2 = trips[i].chi2;
1229 // best_pos = i;
1230 // continue;
1231 // }
1232 // if ((trips[i].chi2 < best_chi2) ||
1233 // (usetrack[trips[best_pos].track] == false)) {
1234 // next_chi2 = best_chi2;
1235 // best_chi2 = trips[i].chi2;
1236 // best_pos = i;
1237 // } else if (trips[i].chi2 < next_chi2) {
1238 // next_chi2 = trips[i].chi2;
1239 // }
1240 // }
1241 // for (unsigned int i = pos; i < next_pos; ++i) {
1242 // if (i != best_pos) {
1243 // usetrack[trips[i].track] = false;
1244 // }
1245 // }
1246 // }
1247 // pos = next_pos;
1248 // cur_h1 = trips[pos].hit1;
1249 // cur_h2 = trips[pos].hit2;
1250 // }
1251 // }
1252 
1253 void sPHENIXTracker::findTracksBySegments(vector<SimpleHit3D>& hits,
1254  vector<SimpleTrack3D>& tracks,
1255  const HelixRange& /*range*/) {
1256 
1257  cout << "fidTracksBySegments " << endl;
1258  vector<TrackSegment>* cur_seg = &segments1;
1259  vector<TrackSegment>* next_seg = &segments2;
1260  unsigned int curseg_size = 0;
1261  unsigned int nextseg_size = 0;
1262 
1263  vector<TrackSegment> complete_segments;
1264 
1265  unsigned int allowed_missing = n_layers - req_layers;
1266 
1267  for (unsigned int l = 0; l < n_layers; ++l) {
1268  layer_sorted[l].clear();
1269  }
1270  for (unsigned int i = 0; i < hits.size(); ++i) {
1271  unsigned int min = (hits[i].get_layer() - allowed_missing);
1272  if (allowed_missing > (unsigned int) hits[i].get_layer()) {
1273  min = 0;
1274  }
1275  for (int l = min; l <= hits[i].get_layer(); l += 1) {
1276  layer_sorted[l].push_back(hits[i]);
1277  }
1278  }
1279  for (unsigned int l = 0; l < n_layers; ++l) {
1280  if (layer_sorted[l].size() == 0) {
1281  return;
1282  }
1283  }
1284 
1285  timeval t1, t2;
1286  double time1 = 0.;
1287  double time2 = 0.;
1288 
1289  gettimeofday(&t1, NULL);
1290 
1291  float cosang_diff = 1. - cosang_cut;
1292  float cosang_diff_inv = 1. / cosang_diff;
1293  float sinang_cut = sqrt(1. - cosang_cut * cosang_cut);
1294  float easy_chi2_cut = ca_chi2_cut;
1295 
1296  vector<float> inv_layer;
1297  inv_layer.assign(n_layers, 1.);
1298  for (unsigned int l = 3; l < n_layers; ++l) {
1299  inv_layer[l] = 1. / (((float)l) - 2.);
1300  }
1301 
1302  unsigned int hit_counter = 0;
1303  float x1_a[4] __attribute__((aligned(16)));
1304  float x2_a[4] __attribute__((aligned(16)));
1305  float x3_a[4] __attribute__((aligned(16)));
1306  float y1_a[4] __attribute__((aligned(16)));
1307  float y2_a[4] __attribute__((aligned(16)));
1308  float y3_a[4] __attribute__((aligned(16)));
1309  float z1_a[4] __attribute__((aligned(16)));
1310  float z2_a[4] __attribute__((aligned(16)));
1311  float z3_a[4] __attribute__((aligned(16)));
1312 
1313  float dx1_a[4] __attribute__((aligned(16)));
1314  float dx2_a[4] __attribute__((aligned(16)));
1315  float dx3_a[4] __attribute__((aligned(16)));
1316  float dy1_a[4] __attribute__((aligned(16)));
1317  float dy2_a[4] __attribute__((aligned(16)));
1318  float dy3_a[4] __attribute__((aligned(16)));
1319  float dz1_a[4] __attribute__((aligned(16)));
1320  float dz2_a[4] __attribute__((aligned(16)));
1321  float dz3_a[4] __attribute__((aligned(16)));
1322 
1323  float kappa_a[4] __attribute__((aligned(16)));
1324  float dkappa_a[4] __attribute__((aligned(16)));
1325 
1326  float ux_mid_a[4] __attribute__((aligned(16)));
1327  float uy_mid_a[4] __attribute__((aligned(16)));
1328  float ux_end_a[4] __attribute__((aligned(16)));
1329  float uy_end_a[4] __attribute__((aligned(16)));
1330 
1331  float dzdl_1_a[4] __attribute__((aligned(16)));
1332  float dzdl_2_a[4] __attribute__((aligned(16)));
1333  float ddzdl_1_a[4] __attribute__((aligned(16)));
1334  float ddzdl_2_a[4] __attribute__((aligned(16)));
1335 
1336  float cur_kappa_a[4] __attribute__((aligned(16)));
1337  float cur_dkappa_a[4] __attribute__((aligned(16)));
1338  float cur_ux_a[4] __attribute__((aligned(16)));
1339  float cur_uy_a[4] __attribute__((aligned(16)));
1340  float cur_chi2_a[4] __attribute__((aligned(16)));
1341  float chi2_a[4] __attribute__((aligned(16)));
1342 
1343  unsigned int hit1[4];
1344  unsigned int hit2[4];
1345  unsigned int hit3[4];
1346 
1347  TrackSegment temp_segment;
1348  temp_segment.hits.assign(n_layers, 0);
1349  // make segments out of first 3 layers
1350  for (unsigned int i = 0, sizei = layer_sorted[0].size(); i < sizei; ++i) {
1351  for (unsigned int j = 0, sizej = layer_sorted[1].size(); j < sizej; ++j) {
1352  for (unsigned int k = 0, sizek = layer_sorted[2].size(); k < sizek; ++k) {
1353  if ((layer_sorted[0][i].get_layer() >= layer_sorted[1][j].get_layer()) ||
1354  (layer_sorted[1][j].get_layer() >= layer_sorted[2][k].get_layer())) {
1355  continue;
1356  }
1357 
1358  x1_a[hit_counter] = layer_sorted[0][i].get_x();
1359  y1_a[hit_counter] = layer_sorted[0][i].get_y();
1360  z1_a[hit_counter] = layer_sorted[0][i].get_z();
1361 
1363 
1364  dx1_a[hit_counter] = 0.5*sqrt(12.0)*sqrt(layer_sorted[0][i].get_size(0,0));
1365  dy1_a[hit_counter] = 0.5*sqrt(12.0)*sqrt(layer_sorted[0][i].get_size(1,1));
1366  dz1_a[hit_counter] = 0.5*sqrt(12.0)*sqrt(layer_sorted[0][i].get_size(2,2));
1367 
1368  x2_a[hit_counter] = layer_sorted[1][j].get_x();
1369  y2_a[hit_counter] = layer_sorted[1][j].get_y();
1370  z2_a[hit_counter] = layer_sorted[1][j].get_z();
1371 
1372  dx2_a[hit_counter] = 0.5*sqrt(12.0)*sqrt(layer_sorted[1][j].get_size(0,0));
1373  dy2_a[hit_counter] = 0.5*sqrt(12.0)*sqrt(layer_sorted[1][j].get_size(1,1));
1374  dz2_a[hit_counter] = 0.5*sqrt(12.0)*sqrt(layer_sorted[1][j].get_size(2,2));
1375 
1376  x3_a[hit_counter] = layer_sorted[2][k].get_x();
1377  y3_a[hit_counter] = layer_sorted[2][k].get_y();
1378  z3_a[hit_counter] = layer_sorted[2][k].get_z();
1379 
1380  dx3_a[hit_counter] = 0.5*sqrt(12.0)*sqrt(layer_sorted[2][k].get_size(0,0));
1381  dy3_a[hit_counter] = 0.5*sqrt(12.0)*sqrt(layer_sorted[2][k].get_size(1,1));
1382  dz3_a[hit_counter] = 0.5*sqrt(12.0)*sqrt(layer_sorted[2][k].get_size(2,2));
1383 
1384  hit1[hit_counter] = i;
1385  hit2[hit_counter] = j;
1386  hit3[hit_counter] = k;
1387 
1388  hit_counter += 1;
1389 
1390  if (hit_counter == 4) {
1391  calculateKappaTangents(x1_a, y1_a, z1_a, x2_a, y2_a, z2_a, x3_a, y3_a,
1392  z3_a, dx1_a, dy1_a, dz1_a, dx2_a, dy2_a, dz2_a,
1393  dx3_a, dy3_a, dz3_a, kappa_a, dkappa_a,
1394  ux_mid_a, uy_mid_a, ux_end_a, uy_end_a,
1395  dzdl_1_a, dzdl_2_a, ddzdl_1_a, ddzdl_2_a);
1396 
1397  for (unsigned int h = 0; h < hit_counter; ++h) {
1398  temp_segment.chi2 =
1399  (dzdl_1_a[h] - dzdl_2_a[h]) /
1400  (ddzdl_1_a[h] + ddzdl_2_a[h] + fabs(dzdl_1_a[h] * sinang_cut));
1401  temp_segment.chi2 *= temp_segment.chi2;
1402  if (temp_segment.chi2 > 2.0) {
1403  continue;
1404  }
1405  temp_segment.ux = ux_end_a[h];
1406  temp_segment.uy = uy_end_a[h];
1407  temp_segment.kappa = kappa_a[h];
1408  if (temp_segment.kappa > top_range.max_k) {
1409  continue;
1410  }
1411  temp_segment.dkappa = dkappa_a[h];
1412  temp_segment.hits[0] = hit1[h];
1413  temp_segment.hits[1] = hit2[h];
1414  temp_segment.hits[2] = hit3[h];
1415  temp_segment.n_hits = 3;
1416  unsigned int outer_layer =
1417  layer_sorted[2][temp_segment.hits[2]].get_layer();
1418  if ((outer_layer - 2) > allowed_missing) {
1419  continue;
1420  }
1421  if ((n_layers - 3) <= allowed_missing) {
1422  complete_segments.push_back(temp_segment);
1423  }
1424  if (next_seg->size() == nextseg_size) {
1425  next_seg->push_back(temp_segment);
1426  nextseg_size += 1;
1427  } else {
1428  (*next_seg)[nextseg_size] = temp_segment;
1429  nextseg_size += 1;
1430  }
1431  }
1432 
1433  hit_counter = 0;
1434  }
1435  }
1436  }
1437  }
1438  if (hit_counter != 0) {
1439  calculateKappaTangents(x1_a, y1_a, z1_a, x2_a, y2_a, z2_a, x3_a, y3_a, z3_a,
1440  dx1_a, dy1_a, dz1_a, dx2_a, dy2_a, dz2_a, dx3_a,
1441  dy3_a, dz3_a, kappa_a, dkappa_a, ux_mid_a, uy_mid_a,
1442  ux_end_a, uy_end_a, dzdl_1_a, dzdl_2_a, ddzdl_1_a,
1443  ddzdl_2_a);
1444 
1445  for (unsigned int h = 0; h < hit_counter; ++h) {
1446  temp_segment.chi2 =
1447  (dzdl_1_a[h] - dzdl_2_a[h]) /
1448  (ddzdl_1_a[h] + ddzdl_2_a[h] + fabs(dzdl_1_a[h] * sinang_cut));
1449  temp_segment.chi2 *= temp_segment.chi2;
1450  if (temp_segment.chi2 > 2.0) {
1451  continue;
1452  }
1453  temp_segment.ux = ux_end_a[h];
1454  temp_segment.uy = uy_end_a[h];
1455  temp_segment.kappa = kappa_a[h];
1456  if (temp_segment.kappa > top_range.max_k) {
1457  continue;
1458  }
1459  temp_segment.dkappa = dkappa_a[h];
1460  temp_segment.hits[0] = hit1[h];
1461  temp_segment.hits[1] = hit2[h];
1462  temp_segment.hits[2] = hit3[h];
1463  temp_segment.n_hits = 3;
1464  unsigned int outer_layer = layer_sorted[2][temp_segment.hits[2]].get_layer();
1465  if ((outer_layer - 2) > allowed_missing) {
1466  continue;
1467  }
1468  if ((n_layers - 3) <= allowed_missing) {
1469  complete_segments.push_back(temp_segment);
1470  }
1471  if (next_seg->size() == nextseg_size) {
1472  next_seg->push_back(temp_segment);
1473  nextseg_size += 1;
1474  } else {
1475  (*next_seg)[nextseg_size] = temp_segment;
1476  nextseg_size += 1;
1477  }
1478  }
1479 
1480  hit_counter = 0;
1481  }
1482 
1483  swap(cur_seg, next_seg);
1484  swap(curseg_size, nextseg_size);
1485 
1486  // add hits to segments layer-by-layer, cutting out bad segments
1487  unsigned int whichseg[4];
1488  for (unsigned int l = 3; l < n_layers; ++l) {
1489  if (l == (n_layers - 1)) {
1490  easy_chi2_cut *= 0.25;
1491  }
1492  nextseg_size = 0;
1493  for (unsigned int i = 0, sizei = curseg_size; i < sizei; ++i) {
1494  for (unsigned int j = 0, sizej = layer_sorted[l].size(); j < sizej; ++j) {
1495  if ((layer_sorted[l - 1][(*cur_seg)[i].hits[l - 1]].get_layer() >=
1496  layer_sorted[l][j].get_layer())) {
1497  continue;
1498  }
1499 
1500  x1_a[hit_counter] = layer_sorted[l - 2][(*cur_seg)[i].hits[l - 2]].get_x();
1501  y1_a[hit_counter] = layer_sorted[l - 2][(*cur_seg)[i].hits[l - 2]].get_y();
1502  z1_a[hit_counter] = layer_sorted[l - 2][(*cur_seg)[i].hits[l - 2]].get_z();
1503  x2_a[hit_counter] = layer_sorted[l - 1][(*cur_seg)[i].hits[l - 1]].get_x();
1504  y2_a[hit_counter] = layer_sorted[l - 1][(*cur_seg)[i].hits[l - 1]].get_y();
1505  z2_a[hit_counter] = layer_sorted[l - 1][(*cur_seg)[i].hits[l - 1]].get_z();
1506  x3_a[hit_counter] = layer_sorted[l][j].get_x();
1507  y3_a[hit_counter] = layer_sorted[l][j].get_y();
1508  z3_a[hit_counter] = layer_sorted[l][j].get_z();
1509 
1510  dx1_a[hit_counter] = 0.5*sqrt(12.0)*sqrt(layer_sorted[l - 2][(*cur_seg)[i].hits[l - 2]].get_size(0,0));
1511  dy1_a[hit_counter] = 0.5*sqrt(12.0)*sqrt(layer_sorted[l - 2][(*cur_seg)[i].hits[l - 2]].get_size(1,1));
1512  dz1_a[hit_counter] = 0.5*sqrt(12.0)*sqrt(layer_sorted[l - 2][(*cur_seg)[i].hits[l - 2]].get_size(2,2));
1513  dx2_a[hit_counter] = 0.5*sqrt(12.0)*sqrt(layer_sorted[l - 1][(*cur_seg)[i].hits[l - 1]].get_size(0,0));
1514  dy2_a[hit_counter] = 0.5*sqrt(12.0)*sqrt(layer_sorted[l - 1][(*cur_seg)[i].hits[l - 1]].get_size(1,1));
1515  dz2_a[hit_counter] = 0.5*sqrt(12.0)*sqrt(layer_sorted[l - 1][(*cur_seg)[i].hits[l - 1]].get_size(2,2));
1516  dx3_a[hit_counter] = 0.5*sqrt(12.0)*sqrt(layer_sorted[l][j].get_size(0,0));
1517  dy3_a[hit_counter] = 0.5*sqrt(12.0)*sqrt(layer_sorted[l][j].get_size(1,1));
1518  dz3_a[hit_counter] = 0.5*sqrt(12.0)*sqrt(layer_sorted[l][j].get_size(2,2));
1519 
1520  cur_kappa_a[hit_counter] = (*cur_seg)[i].kappa;
1521  cur_dkappa_a[hit_counter] = (*cur_seg)[i].dkappa;
1522  cur_ux_a[hit_counter] = (*cur_seg)[i].ux;
1523  cur_uy_a[hit_counter] = (*cur_seg)[i].uy;
1524  cur_chi2_a[hit_counter] = (*cur_seg)[i].chi2;
1525 
1526  whichseg[hit_counter] = i;
1527  hit1[hit_counter] = j;
1528 
1529  hit_counter += 1;
1530  if (hit_counter == 4) {
1532  x1_a, y1_a, z1_a, x2_a, y2_a, z2_a, x3_a, y3_a, z3_a, dx1_a,
1533  dy1_a, dz1_a, dx2_a, dy2_a, dz2_a, dx3_a, dy3_a, dz3_a, kappa_a,
1534  dkappa_a, ux_mid_a, uy_mid_a, ux_end_a, uy_end_a, dzdl_1_a,
1535  dzdl_2_a, ddzdl_1_a, ddzdl_2_a, sinang_cut, cosang_diff_inv,
1536  cur_kappa_a, cur_dkappa_a, cur_ux_a, cur_uy_a, cur_chi2_a,
1537  chi2_a);
1538 
1539  for (unsigned int h = 0; h < hit_counter; ++h) {
1540  if ((chi2_a[h]) * inv_layer[l] < easy_chi2_cut) {
1541  temp_segment.chi2 = chi2_a[h];
1542  temp_segment.ux = ux_end_a[h];
1543  temp_segment.uy = uy_end_a[h];
1544  temp_segment.kappa = kappa_a[h];
1545  if (temp_segment.kappa > top_range.max_k) {
1546  continue;
1547  }
1548  temp_segment.dkappa = dkappa_a[h];
1549  for (unsigned int ll = 0; ll < l; ++ll) {
1550  temp_segment.hits[ll] = (*cur_seg)[whichseg[h]].hits[ll];
1551  }
1552  temp_segment.hits[l] = hit1[h];
1553  unsigned int outer_layer =
1554  layer_sorted[l][temp_segment.hits[l]].get_layer();
1555  temp_segment.n_hits = l + 1;
1556  if ((n_layers - (l + 1)) <= allowed_missing) {
1557  complete_segments.push_back(temp_segment);
1558  }
1559  if ((outer_layer - l) > allowed_missing) {
1560  continue;
1561  }
1562  if (next_seg->size() == nextseg_size) {
1563  next_seg->push_back(temp_segment);
1564  nextseg_size += 1;
1565  } else {
1566  (*next_seg)[nextseg_size] = temp_segment;
1567  nextseg_size += 1;
1568  }
1569  }
1570  }
1571  hit_counter = 0;
1572  }
1573  }
1574  }
1575  if (hit_counter != 0) {
1577  x1_a, y1_a, z1_a, x2_a, y2_a, z2_a, x3_a, y3_a, z3_a, dx1_a, dy1_a,
1578  dz1_a, dx2_a, dy2_a, dz2_a, dx3_a, dy3_a, dz3_a, kappa_a, dkappa_a,
1579  ux_mid_a, uy_mid_a, ux_end_a, uy_end_a, dzdl_1_a, dzdl_2_a, ddzdl_1_a,
1580  ddzdl_2_a, sinang_cut, cosang_diff_inv, cur_kappa_a, cur_dkappa_a,
1581  cur_ux_a, cur_uy_a, cur_chi2_a, chi2_a);
1582 
1583  for (unsigned int h = 0; h < hit_counter; ++h) {
1584  if ((chi2_a[h]) * inv_layer[l] < easy_chi2_cut) {
1585  temp_segment.chi2 = chi2_a[h];
1586  temp_segment.ux = ux_end_a[h];
1587  temp_segment.uy = uy_end_a[h];
1588  temp_segment.kappa = kappa_a[h];
1589  if (temp_segment.kappa > top_range.max_k) {
1590  continue;
1591  }
1592 
1593  temp_segment.dkappa = dkappa_a[h];
1594  for (unsigned int ll = 0; ll < l; ++ll) {
1595  temp_segment.hits[ll] = (*cur_seg)[whichseg[h]].hits[ll];
1596  }
1597  temp_segment.hits[l] = hit1[h];
1598  unsigned int outer_layer =
1599  layer_sorted[l][temp_segment.hits[l]].get_layer();
1600  temp_segment.n_hits = l + 1;
1601  if ((n_layers - (l + 1)) <= allowed_missing) {
1602  complete_segments.push_back(temp_segment);
1603  }
1604  if ((outer_layer - l) > allowed_missing) {
1605  continue;
1606  }
1607  if (next_seg->size() == nextseg_size) {
1608  next_seg->push_back(temp_segment);
1609  nextseg_size += 1;
1610  } else {
1611  (*next_seg)[nextseg_size] = temp_segment;
1612  nextseg_size += 1;
1613  }
1614  }
1615  }
1616  hit_counter = 0;
1617  }
1618  swap(cur_seg, next_seg);
1619  swap(curseg_size, nextseg_size);
1620  }
1621 
1622  for (unsigned int i = 0; i < complete_segments.size(); ++i) {
1623  if (cur_seg->size() == curseg_size) {
1624  cur_seg->push_back(complete_segments[i]);
1625  curseg_size += 1;
1626  } else {
1627  (*cur_seg)[curseg_size] = complete_segments[i];
1628  curseg_size += 1;
1629  }
1630  }
1631 
1632  gettimeofday(&t2, NULL);
1633  time1 = ((double)(t1.tv_sec) + (double)(t1.tv_usec) / 1000000.);
1634  time2 = ((double)(t2.tv_sec) + (double)(t2.tv_usec) / 1000000.);
1635  CAtime += (time2 - time1);
1636  SimpleTrack3D temp_track;
1637  temp_track.hits.assign(n_layers, SimpleHit3D());
1638  vector<SimpleHit3D> temp_hits;
1639  for (unsigned int i = 0, sizei = curseg_size; i < sizei; ++i) {
1640  temp_track.hits.assign((*cur_seg)[i].n_hits, SimpleHit3D());
1641 
1642  temp_comb.assign((*cur_seg)[i].n_hits, 0);
1643  for (unsigned int l = 0; l < (*cur_seg)[i].n_hits; ++l) {
1644  temp_comb[l] = layer_sorted[l][(*cur_seg)[i].hits[l]].get_id();
1645  }
1646  sort(temp_comb.begin(), temp_comb.end());
1647  set<vector<unsigned int> >::iterator it = combos.find(temp_comb);
1648  if (it != combos.end()) {
1649  continue;
1650  }
1651  if (combos.size() > 10000) {
1652  combos.clear();
1653  }
1654  combos.insert(temp_comb);
1655 
1656  for (unsigned int l = 0; l < (*cur_seg)[i].n_hits; ++l) {
1657  temp_track.hits[l] = layer_sorted[l][(*cur_seg)[i].hits[l]];
1658  }
1659 
1660  gettimeofday(&t1, NULL);
1661 
1662  // unsigned int layer_out = temp_track.hits.size()-1;
1663  // unsigned int layer_mid = layer_out/2;
1664  // temp_track_3hits.hits[0] = temp_track.hits[0];
1665  // temp_track_3hits.hits[1] = temp_track.hits[layer_mid];
1666  // temp_track_3hits.hits[2] = temp_track.hits[layer_out];
1667 
1668  float init_chi2 = fitTrack(temp_track);
1669 
1670  if (init_chi2 > fast_chi2_cut_max) {
1671  if (init_chi2 > fast_chi2_cut_par0 +
1672  fast_chi2_cut_par1 / kappaToPt(temp_track.kappa)) {
1673  gettimeofday(&t2, NULL);
1674  time1 = ((double)(t1.tv_sec) + (double)(t1.tv_usec) / 1000000.);
1675  time2 = ((double)(t2.tv_sec) + (double)(t2.tv_usec) / 1000000.);
1676  KALtime += (time2 - time1);
1677  continue;
1678  }
1679  }
1680  HelixKalmanState state;
1681  state.phi = temp_track.phi;
1682  if (state.phi < 0.) {
1683  state.phi += 2. * M_PI;
1684  }
1685  state.d = temp_track.d;
1686  state.kappa = temp_track.kappa;
1687  state.nu = sqrt(state.kappa);
1688  state.z0 = temp_track.z0;
1689  state.dzdl = temp_track.dzdl;
1690  state.C = Matrix<float, 5, 5>::Zero(5, 5);
1691  state.C(0, 0) = pow(0.01, 2.);
1692  state.C(1, 1) = pow(0.01, 2.);
1693  state.C(2, 2) = pow(0.01 * state.nu, 2.);
1694  state.C(3, 3) = pow(0.05, 2.);
1695  state.C(4, 4) = pow(0.05, 2.);
1696  state.chi2 = 0.;
1697  state.position = 0;
1698  state.x_int = 0.;
1699  state.y_int = 0.;
1700  state.z_int = 0.;
1701 
1702  for (unsigned int h = 0; h < temp_track.hits.size(); ++h) {
1703  kalman->addHit(temp_track.hits[h], state);
1704  nfits += 1;
1705  }
1706 
1707  // fudge factor for non-gaussian hit sizes
1708  state.C *= 3.;
1709  state.chi2 *= 6.;
1710 
1711  gettimeofday(&t2, NULL);
1712  time1 = ((double)(t1.tv_sec) + (double)(t1.tv_usec) / 1000000.);
1713  time2 = ((double)(t2.tv_sec) + (double)(t2.tv_usec) / 1000000.);
1714  KALtime += (time2 - time1);
1715 
1716  if (!(temp_track.kappa == temp_track.kappa)) {
1717  continue;
1718  }
1719  if (temp_track.kappa > top_range.max_k) {
1720  continue;
1721  }
1722  if (!(state.chi2 == state.chi2)) {
1723  continue;
1724  }
1725  if (state.chi2 / (2. * ((float)(temp_track.hits.size())) - 5.) > chi2_cut) {
1726  continue;
1727  }
1728 
1729  if (cut_on_dca == true) {
1730  if (fabs(temp_track.d) > dca_cut) {
1731  continue;
1732  }
1733  if (fabs(temp_track.z0) > dca_cut) {
1734  continue;
1735  }
1736  }
1737 
1738  tracks.push_back(temp_track);
1739  track_states.push_back(state);
1740  if ((remove_hits == true) && (state.chi2 < chi2_removal_cut) &&
1741  (temp_track.hits.size() >= n_removal_hits)) {
1742  for (unsigned int i = 0; i < temp_track.hits.size(); ++i) {
1743  (*hit_used)[temp_track.hits[i].get_id()] = true;
1744  }
1745  }
1746  }
1747 }
1748 
1749 
1750 
1751 static inline double sign(double x) {
1752  return ((double)(x > 0.)) - ((double)(x < 0.));
1753 }
1754 
1756  float& x, float& y, float& z) {
1757  float phi = seed.phi;
1758  float d = seed.d;
1759  float k = seed.kappa;
1760  float z0 = seed.z0;
1761  float dzdl = seed.dzdl;
1762 
1763  float hitx = seed.hits.back().get_x();
1764  float hity = seed.hits.back().get_y();
1765 
1766  float rad_det = detector_radii[layer];
1767 
1768  float cosphi = cos(phi);
1769  float sinphi = sin(phi);
1770 
1771  k = fabs(k);
1772 
1773  float kd = (d * k + 1.);
1774  float kcx = kd * cosphi;
1775  float kcy = kd * sinphi;
1776  float kd_inv = 1. / kd;
1777  float R2 = rad_det * rad_det;
1778  float a = 0.5 * (k * R2 + (d * d * k + 2. * d)) * kd_inv;
1779  float tmp1 = a * kd_inv;
1780  float P2x = kcx * tmp1;
1781  float P2y = kcy * tmp1;
1782 
1783  float h = sqrt(R2 - a * a);
1784 
1785  float ux = -kcy * kd_inv;
1786  float uy = kcx * kd_inv;
1787 
1788  float x1 = P2x + ux * h;
1789  float y1 = P2y + uy * h;
1790  float x2 = P2x - ux * h;
1791  float y2 = P2y - uy * h;
1792  float diff1 = (x1 - hitx) * (x1 - hitx) + (y1 - hity) * (y1 - hity);
1793  float diff2 = (x2 - hitx) * (x2 - hitx) + (y2 - hity) * (y2 - hity);
1794  float signk = 0.;
1795  if (diff1 < diff2) {
1796  signk = 1.;
1797  } else {
1798  signk = -1.;
1799  }
1800  x = P2x + signk * ux * h;
1801  y = P2y + signk * uy * h;
1802 
1803  double sign_dzdl = sign(dzdl);
1804 // double onedzdl2_inv = 1. / (1. - dzdl * dzdl);
1805  double startx = d * cosphi;
1806  double starty = d * sinphi;
1807  double D = sqrt((startx - x) * (startx - x) + (starty - y) * (starty - y));
1808 // double D_inv = 1. / D;
1809  double v = 0.5 * k * D;
1810  z = 0.;
1811  if (v > 0.1) {
1812  if (v >= 0.999999) {
1813  v = 0.999999;
1814  }
1815  double s = 2. * asin(v) / k;
1816 // double s_inv = 1. / s;
1817 // double sqrtvv = sqrt(1 - v * v);
1818  double dz = sqrt(s * s * dzdl * dzdl / (1. - dzdl * dzdl));
1819  z = z0 + sign_dzdl * dz;
1820  } else {
1821  double s = 0.;
1822  double temp1 = k * D * 0.5;
1823  temp1 *= temp1;
1824  double temp2 = D * 0.5;
1825  s += 2. * temp2;
1826  temp2 *= temp1;
1827  s += temp2 / 3.;
1828  temp2 *= temp1;
1829  s += (3. / 20.) * temp2;
1830  temp2 *= temp1;
1831  s += (5. / 56.) * temp2;
1832 // double s_inv = 1. / s;
1833  double dz = sqrt(s * s * dzdl * dzdl / (1. - dzdl * dzdl));
1834  z = z0 + sign_dzdl * dz;
1835  }
1836 }
1837 
1838 void sPHENIXTracker::initSplitting(vector<SimpleHit3D>& hits,
1839  unsigned int min_hits,
1840  unsigned int /*max_hits*/) {
1841  initEvent(hits, min_hits);
1842  (*(hits_vec[0])) = hits;
1843  zoomranges.clear();
1844  for (unsigned int z = 0; z <= max_zoom; z++) {
1845  zoomranges.push_back(top_range);
1846  }
1847 }
1848 
1850  vector<SimpleHit3D>& hits, unsigned int /*min_hits*/, unsigned int /*max_hits*/,
1851  vector<SimpleTrack3D>& /*tracks*/) {
1852  unsigned int hits_per_thread = (hits.size() + 2 * nthreads) / nthreads;
1853  unsigned int pos = 0;
1854  while (pos < hits.size()) {
1855  for (unsigned int i = 0; i < nthreads; ++i) {
1856  if (pos >= hits.size()) {
1857  break;
1858  }
1859  for (unsigned int j = 0; j < hits_per_thread; ++j) {
1860  if (pos >= hits.size()) {
1861  break;
1862  }
1863  split_input_hits[i].push_back(hits[pos]);
1864  pos += 1;
1865  }
1866  }
1867  }
1868  for (unsigned int i = 0; i < nthreads; ++i) {
1869  thread_trackers[i]->setTopRange(top_range);
1870  thread_trackers[i]->initSplitting(split_input_hits[i], thread_min_hits,
1871  thread_max_hits);
1872  }
1874  thread_ranges.clear();
1875  thread_hits.clear();
1876 
1877  unsigned int nbins = split_output_hits[0]->size();
1878  for (unsigned int b = 0; b < nbins; ++b) {
1879  thread_ranges.push_back((*(split_ranges[0]))[b]);
1880  thread_hits.push_back(vector<SimpleHit3D>());
1881  for (unsigned int i = 0; i < nthreads; ++i) {
1882  for (unsigned int j = 0; j < (*(split_output_hits[i]))[b].size(); ++j) {
1883  thread_hits.back().push_back((*(split_output_hits[i]))[b][j]);
1884  }
1885  }
1886  }
1887 
1889 }
1890 
1891 void sPHENIXTracker::findHelicesParallel(vector<SimpleHit3D>& hits,
1892  unsigned int min_hits,
1893  unsigned int max_hits,
1894  vector<SimpleTrack3D>& tracks) {
1895  thread_min_hits = min_hits;
1896  thread_max_hits = max_hits;
1897 
1898  for (unsigned int i = 0; i < nthreads; ++i) {
1899  thread_tracks[i].clear();
1900  thread_trackers[i]->clear();
1901  if (cluster_start_bin != 0) {
1902  thread_trackers[i]->setClusterStartBin(cluster_start_bin - 1);
1903  } else {
1904  thread_trackers[i]->setClusterStartBin(0);
1905  }
1906  }
1907 
1908  initSplitting(hits, min_hits, max_hits);
1909 
1910  if (separate_by_helicity == true) {
1911  for (unsigned int i = 0; i < nthreads; ++i) {
1912  thread_trackers[i]->setSeparateByHelicity(true);
1913  thread_trackers[i]->setOnlyOneHelicity(true);
1914  thread_trackers[i]->setHelicity(true);
1915  split_output_hits[i]->clear();
1916  split_input_hits[i].clear();
1917  }
1918  findHelicesParallelOneHelicity(hits, min_hits, max_hits, tracks);
1919 
1920  for (unsigned int i = 0; i < nthreads; ++i) {
1921  thread_trackers[i]->setSeparateByHelicity(true);
1922  thread_trackers[i]->setOnlyOneHelicity(true);
1923  thread_trackers[i]->setHelicity(false);
1924  split_output_hits[i]->clear();
1925  split_input_hits[i].clear();
1926  }
1927  findHelicesParallelOneHelicity(hits, min_hits, max_hits, tracks);
1928  } else {
1929  for (unsigned int i = 0; i < nthreads; ++i) {
1930  thread_trackers[i]->setSeparateByHelicity(false);
1931  thread_trackers[i]->setOnlyOneHelicity(false);
1932  split_output_hits[i]->clear();
1933  split_input_hits[i].clear();
1934  }
1935 
1936  findHelicesParallelOneHelicity(hits, min_hits, max_hits, tracks);
1937  }
1938 
1939  vector<SimpleTrack3D> temp_tracks;
1940  for (unsigned int i = 0; i < nthreads; ++i) {
1941  vector<HelixKalmanState>* states = &(thread_trackers[i]->getKalmanStates());
1942  for (unsigned int j = 0; j < thread_tracks[i].size(); ++j) {
1943  track_states.push_back((*states)[j]);
1944  temp_tracks.push_back(thread_tracks[i][j]);
1945  }
1946  }
1947  finalize(temp_tracks, tracks);
1948 }
1949 
1951  unsigned long int w = (*((unsigned long int*)arg));
1953  *(split_ranges[w]), *(split_output_hits[w]),
1954  0);
1955 }
1956 
1958  unsigned long int w = (*((unsigned long int*)arg));
1959 
1960  for (unsigned int i = w; i < thread_ranges.size(); i += nthreads) {
1961  if (thread_hits[i].size() == 0) {
1962  continue;
1963  }
1964  thread_trackers[w]->setTopRange(thread_ranges[i]);
1965  thread_trackers[w]->findHelices(thread_hits[i], thread_min_hits,
1967  }
1968 }