ECCE @ EIC Software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
VtxTrackFinder.cpp
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file VtxTrackFinder.cpp
1 #include "VtxTrackFinder.h"
2 #include <cmath>
3 #include <iostream>
4 #include <algorithm>
5 #include <Eigen/LU>
6 #include <Eigen/Core>
7 
8 
9 using namespace std;
10 using namespace Eigen;
11 
12 /*
13 bool SimpleHit3D_LessThan(SimpleHit3D i, SimpleHit3D j)
14 {
15  return (i.index<j.index);
16 }
17 
18 bool SimpleTrack3D_Equality(SimpleTrack3D i, SimpleTrack3D j)
19 {
20  if(i.hits.size()!=j.hits.size())
21  return false;
22  for(unsigned int ihit = 0; ihit < i.hits.size(); ihit++)
23  if(i.hits[ihit].index!=j.hits[ihit].index)
24  return false;
25  return true;
26 }
27 
28 bool is_lessthan(vector<SimpleHit3D> i, vector<SimpleHit3D> j, unsigned int ihit)
29 {
30  if(i[ihit].index<j[ihit].index){//if any index for i is < j then it is less than
31  return true;
32  }
33  else if(i[ihit].index==j[ihit].index){//if equal
34  if(ihit < i.size()-1)//and not at the end of the vector
35  return is_lessthan(i,j,ihit+1);//look at the next index
36  else return false;//otherwise at end of vector being equal so return false
37  }
38  else return false;//if index is not < or == then it must be > so return false
39 }
40 
41 bool SimpleTrack3D_LessThan(SimpleTrack3D i , SimpleTrack3D j)
42 {//wrapper for recursive is_lessthan
43 
44 if(i.hits.size()!=j.hits.size()){//if the hits aren't the same size
45  if(i.hits.size()<j.hits.size())//and if no. of hits in i < no. of hits in j
46  return true; //consider i < j
47  else return false; // otherwise i > j
48 }
49 return is_lessthan(i.hits,j.hits,0);//i.size should == j.size so we should be able to compare them index by index
50 
51 }
52 */
53 
54 VtxTrackFinder::VtxTrackFinder(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.)
55 {
56  for(unsigned int i=0;i<detrad.size();++i)
57  {
58  detector_radii.push_back(detrad[i]);
59  }
60 }
61 
62 
63 void VtxTrackFinder::findTracks(vector<SimpleHit3D>& hits, vector<SimpleTrack3D>& tracks)
64 {
65  vector<double> chi2_hits;
66  SimpleTrack3D temp_track;
67  temp_track.hits.resize(4, hits[0]);
68  for(unsigned int i1=0;i1<hits.size();++i1)
69  {
70  temp_track.hits[0] = hits[i1];
71  for(unsigned int i2=(i1+1);i2<hits.size();++i2)
72  {
73  if( (hits[i2].layer == hits[i1].layer)){continue;}
74  temp_track.hits[1] = hits[i2];
75  for(unsigned int i3=(i2+1);i3<hits.size();++i3)
76  {
77  if((hits[i3].layer == hits[i2].layer) || (hits[i3].layer == hits[i1].layer)){continue;}
78  temp_track.hits[2] = hits[i3];
79  for(unsigned int i4=(i3+1);i4<hits.size();++i4)
80  {
81  if( (hits[i4].layer == hits[i3].layer) || (hits[i4].layer == hits[i2].layer) || (hits[i4].layer == hits[i1].layer)){continue;}
82  temp_track.hits[3] = hits[i4];
83 
84  vector<unsigned int> tempcomb;
85  tempcomb.assign(4,0);
86  tempcomb[0] = temp_track.hits[0].index;
87  tempcomb[1] = temp_track.hits[1].index;
88  tempcomb[2] = temp_track.hits[2].index;
89  tempcomb[3] = temp_track.hits[3].index;
90  sort(tempcomb.begin(), tempcomb.end());
91  set<vector<unsigned int> >::iterator it = combos.find(tempcomb);
92  if(it != combos.end()){continue;}
93  combos.insert(tempcomb);
94  double chi2 = fitTrack(temp_track, chi2_hits);
95  if(chi2 > chi2_cut){continue;}
96  tracks.push_back(temp_track);
97  }
98  }
99  }
100  }
101 }
102 
103 
104 void VtxTrackFinder::finalize(vector<SimpleTrack3D>& input, vector<SimpleTrack3D>& output)
105 {
106  unsigned int nt = input.size();
107  for(unsigned int i=0;i<nt;++i)
108  {
109  output.push_back(input[i]);
110  }
111 }
112 
113 
114 double VtxTrackFinder::fitTrack(SimpleTrack3D& track, vector<double>& chi2_hit)
115 {
116  if(using_vertex == true)
117  {
118  track.hits.push_back(SimpleHit3D(0.,0., 0.,0., 0.,0., 0, 0));
119  }
120 
121  chi2_hit.clear();
122  chi2_hit.resize(track.hits.size(), 0.);
123 
124  MatrixXf y = MatrixXf::Zero(track.hits.size(), 1);
125  for(unsigned int i=0;i<track.hits.size();i++)
126  {
127  y(i, 0) = ( pow(track.hits[i].x,2) + pow(track.hits[i].y,2) );
128  if((using_vertex==true ) && (i == (track.hits.size() - 1))){y(i, 0) /= vertex_sigma_xy;}
129  else{y(i, 0) /= layer_xy_resolution[track.hits[i].layer];}
130  }
131 
132  MatrixXf X = MatrixXf::Zero(track.hits.size(), 3);
133  for(unsigned int i=0;i<track.hits.size();i++)
134  {
135  X(i, 0) = track.hits[i].x;
136  X(i, 1) = track.hits[i].y;
137  X(i, 2) = -1.;
138  if((using_vertex==true ) && (i == (track.hits.size() - 1)))
139  {
140  X(i, 0) /= vertex_sigma_xy;
141  X(i, 1) /= vertex_sigma_xy;
142  X(i, 2) /= vertex_sigma_xy;
143  }
144  else
145  {
146  X(i, 0) /= layer_xy_resolution[track.hits[i].layer];
147  X(i, 1) /= layer_xy_resolution[track.hits[i].layer];
148  X(i, 2) /= layer_xy_resolution[track.hits[i].layer];
149  }
150  }
151 
152  MatrixXf Xt = X.transpose();
153 
154  MatrixXf prod = Xt*X;
155  MatrixXf inv = prod.fullPivLu().inverse();
156 
157  MatrixXf beta = inv*Xt*y;
158 
159  float cx = beta(0,0)*0.5;
160  float cy = beta(1,0)*0.5;
161  float r = sqrt(cx*cx + cy*cy - beta(2,0));
162 
163  float phi = atan2(cy, cx);
164  float d = sqrt(cx*cx + cy*cy) - r;
165  float k = 1./r;
166 
167  MatrixXf diff = y - (X*beta);
168  MatrixXf chi2 = (diff.transpose())*diff;
169 
170  float dx = d*cos(phi);
171  float dy = d*sin(phi);
172 
173  MatrixXf y2 = MatrixXf::Zero(track.hits.size(), 1);
174  for(unsigned int i=0;i<track.hits.size();i++)
175  {
176  y2(i,0) = track.hits[i].z;
177  if((using_vertex==true ) && (i == (track.hits.size() - 1))){y2(i, 0) /= vertex_sigma_z;}
178  else{y2(i, 0) /= layer_z_resolution[track.hits[i].layer];}
179  }
180 
181  MatrixXf X2 = MatrixXf::Zero(track.hits.size(), 2);
182  for(unsigned int i=0;i<track.hits.size();i++)
183  {
184  float D = sqrt( pow(dx - track.hits[i].x, 2) + pow(dy - track.hits[i].y,2));
185  float s = 0.0;
186 
187  if(0.5*k*D > 0.1)
188  {
189  float v = 0.5*k*D;
190  if(v >= 0.999999){v = 0.999999;}
191  s = 2.*asin(v)/k;
192  }
193  else
194  {
195  float temp1 = k*D*0.5;temp1*=temp1;
196  float temp2 = D*0.5;
197  s += 2.*temp2;
198  temp2*=temp1;
199  s += temp2/3.;
200  temp2*=temp1;
201  s += (3./20.)*temp2;
202  temp2*=temp1;
203  s += (5./56.)*temp2;
204  }
205 
206  X2(i,0) = s;
207  X2(i,1) = 1.0;
208 
209  if((using_vertex==true ) && (i == (track.hits.size() - 1)))
210  {
211  X2(i, 0) /= vertex_sigma_z;
212  X2(i, 1) /= vertex_sigma_z;
213  }
214  else
215  {
216  X2(i, 0) /= layer_z_resolution[track.hits[i].layer];
217  X2(i, 1) /= layer_z_resolution[track.hits[i].layer];
218  }
219  }
220 
221  MatrixXf Xt2 = X2.transpose();
222  MatrixXf prod2 = Xt2*X2;
223  MatrixXf inv2 = prod2.fullPivLu().inverse();
224  MatrixXf beta2 = inv2*Xt2*y2;
225 
226  MatrixXf diff2 = y2 - (X2*beta2);
227  MatrixXf chi2_z = (diff2.transpose())*diff2;
228 
229  float z0 = beta2(1,0);
230  float dzdl = beta2(0,0)/sqrt(1. + beta2(0,0)*beta2(0,0));
231 
232  track.phi = phi;
233  track.d = d;
234  track.kappa = k;
235  track.dzdl = dzdl;
236  track.z0 = z0;
237 
238  if(track.kappa!=0.)
239  {
240  r=1./track.kappa;
241  }
242  else
243  {
244  r=1.0e10;
245  }
246 
247  cx = (track.d+r)*cos(track.phi);
248  cy = (track.d+r)*sin(track.phi);
249 
250  float chi2_tot = 0.;
251  for(unsigned int h=0;h<track.hits.size();h++)
252  {
253  float dx1 = track.hits[h].x - cx;
254  float dy1 = track.hits[h].y - cy;
255 
256  float dx2 = track.hits[h].x + cx;
257  float dy2 = track.hits[h].y + cy;
258 
259  float xydiff1 = sqrt(dx1*dx1 + dy1*dy1) - r;
260  float xydiff2 = sqrt(dx2*dx2 + dy2*dy2) - r;
261  float xydiff = xydiff2;
262  if(fabs(xydiff1) < fabs(xydiff2)){ xydiff = xydiff1; }
263 
264  float ls_xy = layer_xy_resolution[track.hits[h].layer];
265  if((using_vertex == true) && (h == (track.hits.size() - 1)))
266  {
267  ls_xy = vertex_sigma_xy;
268  }
269 
270  chi2_hit[h] = 0.;
271  chi2_hit[h] += xydiff*xydiff/(ls_xy*ls_xy);
272  chi2_hit[h] += diff2(h,0)*diff2(h,0);
273 
274  chi2_tot += chi2_hit[h];
275  }
276 
277  unsigned int deg_of_freedom = 2*track.hits.size() - 5;
278 
279  if(using_vertex == true)
280  {
281  track.hits.pop_back();
282  chi2_hit.pop_back();
283  }
284 
285  return (chi2_tot)/((double)(deg_of_freedom));
286 }
287 
288