ECCE @ EIC Software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
NHitSeedFinder.cpp
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file NHitSeedFinder.cpp
1 #include "NHitSeedFinder.h"
2 #include <cmath>
3 #include <iostream>
4 #include <algorithm>
5 #include <Eigen/LU>
6 #include <Eigen/Core>
7 
8 using namespace std;
9 using namespace Eigen;
10 
11 NHitSeedFinder::NHitSeedFinder(vector<float>& detrad, unsigned int n_phi,
12  unsigned int n_d, unsigned int n_k,
13  unsigned int n_dzdl, unsigned int n_z0,
14  HelixResolution& min_resolution,
15  HelixResolution& max_resolution, HelixRange& range) :
16  HelixHough(n_phi, n_d, n_k, n_dzdl, n_z0, min_resolution, max_resolution, range),
17  using_vertex(false),
18  vertex_sigma_xy(0.002),
19  vertex_sigma_z(0.005),
20  chi2_cut(3.),
21  seed_mode(4)
22 {
23  for(unsigned int i=0;i<detrad.size();++i)
24  {
25  detector_radii.push_back(detrad[i]);
26  }
27 }
28 
29 void NHitSeedFinder::findTracks(vector<SimpleHit3D>& hits, vector<SimpleTrack3D>& tracks)
30 {
31  if(seed_mode == 4)
32  find4Tracks(hits,tracks);
33  if(seed_mode == 5)
34  find5Tracks(hits,tracks);
35  if(seed_mode == 6)
36  find6Tracks(hits,tracks);
37 }
38 
39 
40 void NHitSeedFinder::find4Tracks(vector<SimpleHit3D>& hits, vector<SimpleTrack3D>& tracks)
41 {
42  // loop over the provided hits
43  for(unsigned int i=0;i<hits.size();i++)
44  {
45  // if this index is not in our map
46  if( phis.find(hits[i].index) == phis.end() )
47  {
48  // calculate the phi now and insert into the map
49  phis.insert(make_pair(hits[i].index,atan2(hits[i].y,hits[i].x)));
50  }
51  }
52 
53  vector<double> chi2_hits;
54  SimpleTrack3D temp_track;
55  temp_track.hits.resize(4, hits[0]);
56 
57  for(unsigned int i1=0;i1<hits.size();++i1)
58  {
59  temp_track.hits[0] = hits[i1];
60 
61  for(unsigned int i2=(i1+1);i2<hits.size();++i2)
62  {
63  if( hits[i2].layer == hits[i1].layer ) continue;
64  if( fabs( phis[hits[i2].index] - phis[hits[i1].index] ) > M_PI_2 ) continue;
65 
66  temp_track.hits[1] = hits[i2];
67 
68  for(unsigned int i3=(i2+1);i3<hits.size();++i3)
69  {
70  if( hits[i3].layer == hits[i2].layer ) continue;
71  if( hits[i3].layer == hits[i1].layer ) continue;
72  if( fabs( phis[hits[i3].index] - phis[hits[i2].index] ) > M_PI_2 ) continue;
73  if( fabs( phis[hits[i3].index] - phis[hits[i1].index] ) > M_PI_2 ) continue;
74 
75  temp_track.hits[2] = hits[i3];
76 
77  for(unsigned int i4=(i3+1);i4<hits.size();++i4)
78  {
79  if( hits[i4].layer == hits[i3].layer ) continue;
80  if( hits[i4].layer == hits[i2].layer ) continue;
81  if( hits[i4].layer == hits[i1].layer ) continue;
82  if( fabs( phis[hits[i4].index] - phis[hits[i3].index] ) > M_PI_2 ) continue;
83  if( fabs( phis[hits[i4].index] - phis[hits[i2].index] ) > M_PI_2 ) continue;
84  if( fabs( phis[hits[i4].index] - phis[hits[i1].index] ) > M_PI_2 ) continue;
85 
86  temp_track.hits[3] = hits[i4];
87 
88  // build track
89  vector<unsigned int> tempcomb;
90  tempcomb.assign(4,0);
91  tempcomb[0] = temp_track.hits[0].index;
92  tempcomb[1] = temp_track.hits[1].index;
93  tempcomb[2] = temp_track.hits[2].index;
94  tempcomb[3] = temp_track.hits[3].index;
95 
96  // check this combo has already been seen
97  sort(tempcomb.begin(), tempcomb.end());
98  set<vector<unsigned int> >::iterator it = combos.find(tempcomb);
99  if(it != combos.end()) continue;
100 
101  // add to list of checked combinations
102  combos.insert(tempcomb);
103 
104  // fit track
105  double chi2 = fitTrack(temp_track, chi2_hits);
106  if(chi2 > chi2_cut) continue;
107 
108  // add track
109  tracks.push_back(temp_track);
110  }
111  }
112  }
113  }
114 }
115 
116 
117 void NHitSeedFinder::find5Tracks(vector<SimpleHit3D>& hits, vector<SimpleTrack3D>& tracks)
118 {//cout<<"find6Tracks Entered with "<<hits.size()<<" hits"<<endl;
119  // loop over the provided hits
120  for(unsigned int i=0;i<hits.size();i++)
121  {
122  // if this index is not in our map
123  if( phis.find(hits[i].index) == phis.end() )
124  {
125  // calculate the phi now and insert into the map
126  phis.insert(make_pair(hits[i].index,atan2(hits[i].y,hits[i].x)));
127  }
128  }
129 
130  vector<double> chi2_hits;
131  SimpleTrack3D temp_track;
132  temp_track.hits.resize(5, hits[0]);
133 
134  for(unsigned int i1=0;i1<hits.size();++i1)
135  {
136  temp_track.hits[0] = hits[i1];
137 
138  for(unsigned int i2=(i1+1);i2<hits.size();++i2)
139  {
140  if( hits[i2].layer == hits[i1].layer ) continue;
141  /* if( fabs( phis[hits[i3].index] - phis[hits[i2].index] ) > M_PI_2 ) continue;
142  if( fabs( phis[hits[i3].index] - phis[hits[i1].index] ) > M_PI_2 ) continue;*/
143 
144  temp_track.hits[1] = hits[i2];
145 
146  for(unsigned int i3=(i2+1);i3<hits.size();++i3)
147  {
148  if( hits[i3].layer == hits[i2].layer ) continue;
149  if( hits[i3].layer == hits[i1].layer ) continue;
150  /* if( fabs( phis[hits[i4].index] - phis[hits[i3].index] ) > M_PI_2 ) continue;
151  if( fabs( phis[hits[i4].index] - phis[hits[i2].index] ) > M_PI_2 ) continue;
152  if( fabs( phis[hits[i4].index] - phis[hits[i1].index] ) > M_PI_2 ) continue;*/
153 
154  temp_track.hits[2] = hits[i3];
155 
156  for(unsigned int i4=(i3+1);i4<hits.size();++i4)
157  {
158  if( hits[i4].layer == hits[i3].layer ) continue;
159  if( hits[i4].layer == hits[i2].layer ) continue;
160  if( hits[i4].layer == hits[i1].layer ) continue;
161  /* if( fabs( phis[hits[i5].index] - phis[hits[i4].index] ) > M_PI_2 ) continue;
162  if( fabs( phis[hits[i5].index] - phis[hits[i3].index] ) > M_PI_2 ) continue;
163  if( fabs( phis[hits[i5].index] - phis[hits[i2].index] ) > M_PI_2 ) continue;
164  if( fabs( phis[hits[i5].index] - phis[hits[i1].index] ) > M_PI_2 ) continue;*/
165 
166  temp_track.hits[3] = hits[i4];
167  for(unsigned int i5=(i4+1);i5<hits.size();++i5)
168  {
169  if( hits[i5].layer == hits[i4].layer ) continue;
170  if( hits[i5].layer == hits[i3].layer ) continue;
171  if( hits[i5].layer == hits[i2].layer ) continue;
172  if( hits[i5].layer == hits[i1].layer ) continue;
173  /* if( fabs( phis[hits[i5].index] - phis[hits[i4].index] ) > M_PI_2 ) continue; if( fabs( phis[hits[i5].index] - phis[hits[i3].index] ) > M_PI_2 ) continue; if( fabs( phis[hits[i5].index] - phis[hits[i2].index] ) > M_PI_2 ) continue; if( fabs( phis[hits[i5].index] - phis[hits[i1].index] ) > M_PI_2 ) continue;*/
174  temp_track.hits[4]= hits[i5];
175  // build track
176 
177  vector<unsigned int> tempcomb;
178  tempcomb.assign(5,0);
179  tempcomb[0] = temp_track.hits[0].index;
180  tempcomb[1] = temp_track.hits[1].index;
181  tempcomb[2] = temp_track.hits[2].index;
182  tempcomb[3] = temp_track.hits[3].index;
183  tempcomb[4] = temp_track.hits[4].index;
184 
185  // check this combo has already been seen
186 
187  sort(tempcomb.begin(), tempcomb.end());
188  set<vector<unsigned int> >::iterator it = combos.find(tempcomb);
189  if(it != combos.end()) continue;
190 
191  // add to list of checked combinations
192 
193  combos.insert(tempcomb);
194 
195  // fit track
196 
197  double chi2 = fitTrack(temp_track, chi2_hits);
198  // cout <<"chi2 of seed track: "<<chi2<<endl;
199  if(chi2 > chi2_cut) continue;
200 
201  // add track
202  tracks.push_back(temp_track);
203  }
204  }
205  }
206  }
207  }
208 }
209 
210 
211 
212 void NHitSeedFinder::find6Tracks(vector<SimpleHit3D>& hits, vector<SimpleTrack3D>& tracks)
213 {
214 // cout << "find6Tracks Entered with " << hits.size()
215 // << " hits, ncombos: " << combos.size() << endl;
216 
217  // loop over the provided hits
218  for(unsigned int i=0;i<hits.size();i++)
219  {
220  // if this index is not in our map
221  if( phis.find(hits[i].index) == phis.end() )
222  {
223  // calculate the phi now and insert into the map
224  phis.insert(make_pair(hits[i].index,atan2(hits[i].y,hits[i].x)));
225  }
226  }
227 
228  vector<double> chi2_hits;
229  SimpleTrack3D temp_track;
230  temp_track.hits.resize(6, hits[0]);
231 
232  for(unsigned int i1=0;i1<hits.size();++i1)
233  {
234  temp_track.hits[0] = hits[i1];
235 
236  for(unsigned int i2=(i1+1);i2<hits.size();++i2)
237  {
238  if( hits[i2].layer == hits[i1].layer ) continue;
239 
240  //if( fabs( phis[hits[i2].index] - phis[hits[i1].index] ) > M_PI_2 ) continue;
241 
242  temp_track.hits[1] = hits[i2];
243 
244  for(unsigned int i3=(i2+1);i3<hits.size();++i3)
245  {
246  if( hits[i3].layer == hits[i2].layer ) continue;
247  if( hits[i3].layer == hits[i1].layer ) continue;
248  /* if( fabs( phis[hits[i3].index] - phis[hits[i2].index] ) > M_PI_2 ) continue;
249  if( fabs( phis[hits[i3].index] - phis[hits[i1].index] ) > M_PI_2 ) continue;*/
250 
251  temp_track.hits[2] = hits[i3];
252 
253  for(unsigned int i4=(i3+1);i4<hits.size();++i4)
254  {
255  if( hits[i4].layer == hits[i3].layer ) continue;
256  if( hits[i4].layer == hits[i2].layer ) continue;
257  if( hits[i4].layer == hits[i1].layer ) continue;
258  /* if( fabs( phis[hits[i4].index] - phis[hits[i3].index] ) > M_PI_2 ) continue;
259  if( fabs( phis[hits[i4].index] - phis[hits[i2].index] ) > M_PI_2 ) continue;
260  if( fabs( phis[hits[i4].index] - phis[hits[i1].index] ) > M_PI_2 ) continue;*/
261 
262  temp_track.hits[3] = hits[i4];
263 
264  for(unsigned int i5=(i4+1);i5<hits.size();++i5)
265  {
266  if( hits[i5].layer == hits[i4].layer ) continue;
267  if( hits[i5].layer == hits[i3].layer ) continue;
268  if( hits[i5].layer == hits[i2].layer ) continue;
269  if( hits[i5].layer == hits[i1].layer ) continue;
270  /* if( fabs( phis[hits[i5].index] - phis[hits[i4].index] ) > M_PI_2 ) continue;
271  if( fabs( phis[hits[i5].index] - phis[hits[i3].index] ) > M_PI_2 ) continue;
272  if( fabs( phis[hits[i5].index] - phis[hits[i2].index] ) > M_PI_2 ) continue;
273  if( fabs( phis[hits[i5].index] - phis[hits[i1].index] ) > M_PI_2 ) continue;*/
274 
275  temp_track.hits[4] = hits[i5];
276 
277  for(unsigned int i6=(i5+1);i6<hits.size();++i6)
278  {
279  if( hits[i6].layer == hits[i5].layer ) continue;
280  if( hits[i6].layer == hits[i4].layer ) continue;
281  if( hits[i6].layer == hits[i3].layer ) continue;
282  if( hits[i6].layer == hits[i2].layer ) continue;
283  if( hits[i6].layer == hits[i1].layer ) continue;
284  /* if( fabs( phis[hits[i6].index] - phis[hits[i5].index] ) > M_PI_2 ) continue;
285  if( fabs( phis[hits[i6].index] - phis[hits[i4].index] ) > M_PI_2 ) continue;
286  if( fabs( phis[hits[i6].index] - phis[hits[i3].index] ) > M_PI_2 ) continue;
287  if( fabs( phis[hits[i6].index] - phis[hits[i2].index] ) > M_PI_2 ) continue;
288  if( fabs( phis[hits[i6].index] - phis[hits[i1].index] ) > M_PI_2 ) continue;*/
289 
290  temp_track.hits[5] = hits[i6];
291 
292  // build track
293  vector<unsigned int> tempcomb;
294  tempcomb.assign(6,0);
295  tempcomb[0] = temp_track.hits[0].index;
296  tempcomb[1] = temp_track.hits[1].index;
297  tempcomb[2] = temp_track.hits[2].index;
298  tempcomb[3] = temp_track.hits[3].index;
299  tempcomb[4] = temp_track.hits[4].index;
300  tempcomb[5] = temp_track.hits[5].index;
301 
302  // check this combo has already been seen
303  sort(tempcomb.begin(), tempcomb.end());
304  set<vector<unsigned int> >::iterator it = combos.find(tempcomb);
305  if(it != combos.end()) {
306  //cout << "combo already used!" << endl;
307  continue;
308  }
309 
310  // add to list of checked combinations
311  combos.insert(tempcomb);
312 
313  // fit track
314  double chi2 = fitTrack(temp_track, chi2_hits);
315  //cout <<"chi2 of seed track: "<<chi2<<endl;
316  if(chi2 > chi2_cut) continue;
317 
318 
319  // add track
320  tracks.push_back(temp_track);
321  }
322  }
323  }
324  }
325  }
326  }
327 }
328 
329 
330 void NHitSeedFinder::finalize(vector<SimpleTrack3D>& input, vector<SimpleTrack3D>& output)
331 {
332  unsigned int nt = input.size();
333  for(unsigned int i=0;i<nt;++i)
334  {
335  output.push_back(input[i]);
336  }
337 }
338 
339 
340 double NHitSeedFinder::fitTrack(SimpleTrack3D& track, vector<double>& chi2_hit)
341 {
342  if(using_vertex == true)
343  {
344  track.hits.push_back(SimpleHit3D(0.,0., 0.,0., 0.,0., 0, 0));
345  }
346 
347  chi2_hit.clear();
348  chi2_hit.resize(track.hits.size(), 0.);
349 
350  MatrixXf y = MatrixXf::Zero(track.hits.size(), 1);
351  for(unsigned int i=0;i<track.hits.size();i++)
352  {
353  y(i, 0) = ( pow(track.hits[i].x,2) + pow(track.hits[i].y,2) );
354  if((using_vertex==true ) && (i == (track.hits.size() - 1))){y(i, 0) /= vertex_sigma_xy;}
355  else{y(i, 0) /= layer_xy_resolution[track.hits[i].layer];}
356  }
357 
358  MatrixXf X = MatrixXf::Zero(track.hits.size(), 3);
359  for(unsigned int i=0;i<track.hits.size();i++)
360  {
361  X(i, 0) = track.hits[i].x;
362  X(i, 1) = track.hits[i].y;
363  X(i, 2) = -1.;
364  if((using_vertex==true ) && (i == (track.hits.size() - 1)))
365  {
366  X(i, 0) /= vertex_sigma_xy;
367  X(i, 1) /= vertex_sigma_xy;
368  X(i, 2) /= vertex_sigma_xy;
369  }
370  else
371  {
372  X(i, 0) /= layer_xy_resolution[track.hits[i].layer];
373  X(i, 1) /= layer_xy_resolution[track.hits[i].layer];
374  X(i, 2) /= layer_xy_resolution[track.hits[i].layer];
375  }
376  }
377 
378  MatrixXf Xt = X.transpose();
379 
380  MatrixXf prod = Xt*X;
381  MatrixXf inv = prod.fullPivLu().inverse();
382 
383  MatrixXf beta = inv*Xt*y;
384 
385  float cx = beta(0,0)*0.5;
386  float cy = beta(1,0)*0.5;
387  float r = sqrt(cx*cx + cy*cy - beta(2,0));
388 
389  float phi = atan2(cy, cx);
390  float d = sqrt(cx*cx + cy*cy) - r;
391  float k = 1./r;
392 
393  MatrixXf diff = y - (X*beta);
394  MatrixXf chi2 = (diff.transpose())*diff;
395 
396  float dx = d*cos(phi);
397  float dy = d*sin(phi);
398 
399  MatrixXf y2 = MatrixXf::Zero(track.hits.size(), 1);
400  for(unsigned int i=0;i<track.hits.size();i++)
401  {
402  y2(i,0) = track.hits[i].z;
403  if((using_vertex==true ) && (i == (track.hits.size() - 1))){y2(i, 0) /= vertex_sigma_z;}
404  else{y2(i, 0) /= layer_z_resolution[track.hits[i].layer];}
405  }
406 
407  MatrixXf X2 = MatrixXf::Zero(track.hits.size(), 2);
408  for(unsigned int i=0;i<track.hits.size();i++)
409  {
410  float D = sqrt( pow(dx - track.hits[i].x, 2) + pow(dy - track.hits[i].y,2));
411  float s = 0.0;
412 
413  if(0.5*k*D > 0.1)
414  {
415  float v = 0.5*k*D;
416  if(v >= 0.999999){v = 0.999999;}
417  s = 2.*asin(v)/k;
418  }
419  else
420  {
421  float temp1 = k*D*0.5;temp1*=temp1;
422  float temp2 = D*0.5;
423  s += 2.*temp2;
424  temp2*=temp1;
425  s += temp2/3.;
426  temp2*=temp1;
427  s += (3./20.)*temp2;
428  temp2*=temp1;
429  s += (5./56.)*temp2;
430  }
431 
432  X2(i,0) = s;
433  X2(i,1) = 1.0;
434 
435  if((using_vertex==true ) && (i == (track.hits.size() - 1)))
436  {
437  X2(i, 0) /= vertex_sigma_z;
438  X2(i, 1) /= vertex_sigma_z;
439  }
440  else
441  {
442  X2(i, 0) /= layer_z_resolution[track.hits[i].layer];
443  X2(i, 1) /= layer_z_resolution[track.hits[i].layer];
444  }
445  }
446 
447  MatrixXf Xt2 = X2.transpose();
448  MatrixXf prod2 = Xt2*X2;
449  MatrixXf inv2 = prod2.fullPivLu().inverse();
450  MatrixXf beta2 = inv2*Xt2*y2;
451 
452  MatrixXf diff2 = y2 - (X2*beta2);
453  MatrixXf chi2_z = (diff2.transpose())*diff2;
454 
455  float z0 = beta2(1,0);
456  float dzdl = beta2(0,0)/sqrt(1. + beta2(0,0)*beta2(0,0));
457 
458  track.phi = phi;
459  track.d = d;
460  track.kappa = k;
461  track.dzdl = dzdl;
462  track.z0 = z0;
463 
464  if(track.kappa!=0.)
465  {
466  r=1./track.kappa;
467  }
468  else
469  {
470  r=1.0e10;
471  }
472 
473  cx = (track.d+r)*cos(track.phi);
474  cy = (track.d+r)*sin(track.phi);
475 
476  float chi2_tot = 0.;
477  for(unsigned int h=0;h<track.hits.size();h++)
478  {
479  float dx1 = track.hits[h].x - cx;
480  float dy1 = track.hits[h].y - cy;
481 
482  float dx2 = track.hits[h].x + cx;
483  float dy2 = track.hits[h].y + cy;
484 
485  float xydiff1 = sqrt(dx1*dx1 + dy1*dy1) - r;
486  float xydiff2 = sqrt(dx2*dx2 + dy2*dy2) - r;
487  float xydiff = xydiff2;
488  if(fabs(xydiff1) < fabs(xydiff2)){ xydiff = xydiff1; }
489 
490  float ls_xy = layer_xy_resolution[track.hits[h].layer];
491  if((using_vertex == true) && (h == (track.hits.size() - 1)))
492  {
493  ls_xy = vertex_sigma_xy;
494  }
495 
496  chi2_hit[h] = 0.;
497  chi2_hit[h] += xydiff*xydiff/(ls_xy*ls_xy);
498  chi2_hit[h] += diff2(h,0)*diff2(h,0);
499 
500  chi2_tot += chi2_hit[h];
501  }
502 
503  unsigned int deg_of_freedom = 2*track.hits.size() - 5;
504 
505  if(using_vertex == true)
506  {
507  track.hits.pop_back();
508  chi2_hit.pop_back();
509  }
510 
511  return (chi2_tot)/((double)(deg_of_freedom));
512 }
513 
514