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